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FINITE STATE MODETTNG OF AFR0P:TASTTC SYSTEMS 


Ran Jan Vepa 
Stanford University 

Sl^MMy\RY 

A general theory of finite state modeling of aerodynamic loads on 
thin airfoils and lifting siir fares performing completely arbitrary, small, 
time-dependent motions in an airstream is . systemati cal ly developed and 
presented. In particular, the nature of the behavior of the unsteady 
airloads in tiie frequency domain is explained. This scheme employs as 
raw materials any of the unsteady linearized theories that have been 
mechanized for simple harmonic oscillations. Each desired aerodynamic 
transfer function is approximated by means of an appropriate Padd 
approximant, that is, a rational function of finite degree polynomials 
in the Laplace transform variable. 

The modeling technique is applied to several two-dimensional and 
three-dimensional airfoils. Circular, elliptic, rectangular and tapered 
planforms are considered as examples. Identical functions are also ob- 
tained for control surfaces for two- and three-dimensional airfoils. 

IN TRODUCTION 

In the last decade rripic] advances have taken place in the area of 
automatic control of practical tongineering systems. The vast tech- 
nological developments in autopilot design and in the design of aircraft 
take-off and landing systems has led to the possibility of applying 
tliis technology to control the vibration mcnies of <iircraft wing structures 
and the elimination of aercu-lastic instabilities in the flight envelope 


of the aircraft. Although the mathematical theory of distributed parameter 
systems has made rapid advances recently, it seems more expedient, from a 
practical point of view, to approximate aeroelastic systems, mainly aircraft 
wings and control surfaces, by finite state models. The techniques of 
approximating aircraft wing structures by finite state models which make use 
of finite elements and other structural idealizations are well known. No 
systematic techniques exist, however, for approximating the aerodynamic loads 
on these structures by compatible finite state models for aeroelastic purposes. 
Thus a systematic theory for approximating aerodynamic loads on aircraft wings 
by finite state models used along with well-known techniques of structural 
idealization could be tremendously useful not only for understanding aero- 
elastic Instabilities but also in the development of control systems for 
suppressing aeroelastic instabilities. Such theories can also prove helpful 
in the minimum weight design of aircraft structures. 

This paper is concerned with the finite state modeling of aeroelastic 
systems. The well-known theories of modeling of aircraft wing structures are 
briefly presented. A general theory is then developed for the modeling of 
imsteady aerodynamic loads on wings and airfoils. These aerodynamic models 
may be used in conjunction with structural models for aeroelastic purposes. 

General Description of the Problem 

The general techniques of calculating unsteady aerodynamic loads for 
simple harmonically oscillating airfoils and lifting surfaces are out- 
lined in [^l] , [2J and [[3]. With these techniques it is possible to 
calculate the unsteady aerodynamic loads for different modes of oscilla- 
tion at a given frequency of oscillation. Little is known about the 
analytical behavior of these loads in the frequency domain. Thus it is 
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first essential to identify the behavior in the frequency domain of the 
aerodynamic loads on airfoils and lifting surfaces. 

The next step is to approximate this behavior in a manner that will 
permit the construction of aerodynamic models, which may be used along 
with structural models for aeroelastic purposes. Also, for the work 
described in this paper existing techniques of calculating aerodynamic 
loads for airfoils and lifting surfaces were utilized whenever possible. 

Review of Pertinent Literature 

The theory of finite state modeling of structures for dynamic 
analysis is well known. The various methods of weighted residuals W. 
finite element techniques [s] and variational techniques [^] have all 
proved extremely useful for analytical modeling purposes. Recently 
system identification techniques [b] have been formulated for approxi- 
mating structures by finite state models from experimental data. 

The author is not aware of any systematic modeling procedures in 
unsteady aerod 3 mamlcs . However, there were several related developments 
in the past fifty years. In 1925, Wagner [l,3] first studied the 
growth of lift on a two dimensional airfoil in incompressible flow due 
to an impulsive change in the vertical velocity of the airfoil. Garrick 
[7] , later showed the relationship between Wagner's solution and Theo- 
dorsen's solution [sj for the lift on an oscillating airfoil. Sears [9] 
showed the relationship between solutions for a sharp edge gust and a 
sinusodial gust. R. T. Jones [10] first considered the aerodynamic forces 
on finite wings of elliptic planform in non-uniform motion in incompressible 
flow. W. P. Jones [llj calculated the lift on rectangular and tapered 
wings for impulsive motion in incompressible flow. Lomax, et al., [l2] 
solved the problem of obtaining the lift and moment for impulsive motion, 
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exactly, at the starting instant tisin:^ i i i >)\v with with supersonic 

steady flow over three dimensional 1 f tt iuu u ■ w; and J)rischler 

[ 13 ] and Drischler [14] used the i-ouri r r r tntis 1 - ^ 1 a t J onsli ip for oscilla- 
tory and impulsive motions developed h. ;;:,e ol tained approximations 

to the indicial lift and moment in pJunpiuu and ri’ jhiiig In ('onip ressible 
flow, over two-dimensional airfoils* dims : . o considered tlie transient 
loading on wide delta wings and rectangUahir vmt'^v o supersonic speeds. 

Recently Djojodihardjo and WidnaM ;17j =uo ii< 0 ; [l^j developed a 
numerical technique for arbitrarily moviue air^d ■ . iu(i lifting surfaces 

based on Green’s representation theorem in I theory. 

While the theory to be presented in rles. -.p r r has some common features 
with each one of the above techniques, one 1^0 pet from them any 

insight into the general behavior of the ae 1 od\ rnir: ^ ^ Icnids and pressure 

distributions in the frequency domain. Tb^ d j t -ut' ces in tiie behavior of 

incompressible and compressible flows have !i(U- h, » li l uJIs^ explained. 

Several techniques have been p>rescaite(i hi r r ,, ra-:t !m>r the calculation 
of aerodynamic loads for simple hannoai' il iv ^ . ; o up airfoils and 

lifting surfaces. Tliese techniques hav^' t)een m v i-w-a bv Ashley, V/idnall 
and Landahl [19] , Landahl and Stark - . r dil j and Ashley ^22]. 

Recent methods now widely in use are tiu ke^-n. : ra t I(m techniques of 
Cunningham [ 2 ^, Rowe [24J and Laschka jr'm, - e tir.. doublet-lattice 
technique of Albano and Redden [26J . iiie . . is \:o he an extremely 

useful and practical tool for aeruelasti ;nn;u s; = . Seme of these 
methods and their extensions are l)riei j v to no m u hi Ap>pendix C. 

In this paper the theory of finite c s Mig of structures is 
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first outlined. Then a systematic procedure is presented for the modeling 
of unsteady aerodynamic loads. Some contributions are also made towards 
extending and modifying methods for calculating aerodynamic loads for 
oscillating airfoils and lifting surfaces. 

Summary of Contributions 

(1) A general theory is presented to describe the behavior of 
the aerodynamic pressure and loads in the frequency domain. 

(2) A numerical technique has been developed for modeling the 
aerodynamic loads in the frequency domain. 

(3) A numerical technique is developed for calculating the aero- 
dynamic loads on two dimensional airfoils with trailing and leading 
edge flaps, for a fairly large frequency bandwidth. For three 
dimensional lifting surfaces the bandwidth, for which the doublet lattice 
method leads to reasonable results, is improved. 

(4) Indlclal functions have been obtained for various types of 
Impulsive motions for wings of different planform. 
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SYMBOLS 


distance of elastic axis from the center of gravity. 


reference semi -chord of wing. 


speed of soiond in the free streani. 


c. 

1 


distance of flap hingeline from the center of gravity, 


c(k) Theodorsen ^ s function . 


non-dimensional pressure coefficient ■: Ap/q 


C (x, I) bending influence function. 

ZZt 


C_„(x,5) torsional influence fixnction. 

C7C7 


El 


flexural rigidity. 


GJ torsional rigidity. 


h(x,t) vertical displacement of vri.ng elastic axis. 


Hankel functions. 


unit matrix. 


moment of inertia of wing. 


P. 

1 


moments of inertia of flaps . 


!.(•) modified Bessel functions of the first kind. 
1 ' ' 


reduced frequency ^ . 



wing root stiffness in shear 



wing root stiffness 

in 

lateral rotation 


wing root stiffness 

in 

torsion. 

K(t) 

Wagner's indlcial function. 


stiffness matrices. 




K^Cx^s) Kernel function for two dimensional aerofoils. 

K2q(X^s) Kernel function for two dimensional aerofoils in incompress- 
ible flow. 

K3 (X,Y,s)Kernel function for three dimensional lifting surfaces. 

K Control surface stiffness constant, 

pi 

defined in figure 1- 
semi-span of wing, 
lift in the plunging mode . 
m(x) mass per unit length. 

Mq,M^ Mass matricies 

M Mach number 

M aerodynamic pitching moment. 

U 

Mq aerodynamic flap moment . 

^i 

NB number of collocation points . 

NF number of flaps 

p Laplace transform variable. 

p(x,s), pressure distribution. 
p(x,y,s) ^ 
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Ap(x,s)^ 
Ap(x ,y,s) 

<1 



(s,M) 

Q(s,M) 


pressure jump across the wing* 

dynamic pressure ~ \f2 pV • 

nodal displacement vector. 

modal displacement vector. 

Laplace transform of q^, q^. 
th 

i, j element of Q(s,M)- 
aerodynamic load matrix- 


s = 


pb 

V 


S,(x) 


“^cpi 


U 


U(n,a,t) 


V 

w(x,y,t) 


unbalance of entire wing, 
th 

imbalance of i^ flap, 
time 

th 

control torque on i^“ flap, 
gust velocity. 

Kummer functions . 


free stream velocity, 
downwash velocity. 
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w(x,y ,p) 


Laplace transform of w(x, y, t) 


x.,y ,z 


t’8,rtesian coordinates. 


X = X - § 


Y y - 'll 


Z = z - C 




til 

displacement of the j mode 


^(s) - s cp(s)‘ 


Vt Vt 

T r- T” 

b g 

w 


angle of attack 


/T 


■ M ■ for M < 1 




1 for M > 1 . 


Pi(t) 


flap rotation of i^^ flap 


0(x,t) 


torsional rotation 


cp(s) 


Laplace transform of Wagner s function- 


c 


Cartesian coordinate system 


FINITE STATE MODELING OF STRUCTURES 


Finite Element Models 

The finite element method is often regarded as generating a dis- 
crete model of a physical system. The method may also be viewed as a 
variational technique, where an attempt is made to minimize the demand 
for ingenuity in the construction of trial functions. This discrete 
matrix technique for the formulation and solution of linear dynamic 
problems in engineering mechanics has been widely discussed in the 
literature [5l* This approach has proved useful in obtaining approxi- 
mate analyses of complex structural configurations that are difficult 
to handle by exact mathematical formulations . The bookkeeping required 
for solving the large number of linear simultaneous equations involved 
is readily handled by matrix algebra techniques and the resulting 
analytical formulations are tremendously simplified. Thus concurrent 
development or revision of different sections of large digital computer 
programs to perform the analysis is feasible. This feature has led to 
its use and acceptance as a basic tool for dynamic structural response 
calculations . The aim of this section is to briefly describe how the 
technique can be useful for finite state modeling of wing structures. 

The efficient utilization of the discrete element influence - 
coefficient approach in a digital computer would provide for con- 
struction of a stiffness matrix and a mass matrix for the entire wing 
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structure by simple superposition of each matrix, from corresponding 
matricies for the discrete elements that model the wing in detail . A £tiffness 
and mass matrix generator program (SAMGEN) was written for this pur- 
pose . Several different types of elements were incorporated in the 
program. This program was found to be extremely useful in structural 

optimization work [27] • 


Idealization as a Beam-Rod 

The swept wing structure shown in Figure 1 is considered. The 
wing is connected structurally to the fuselage , so that the root is not 
completely rigid. Further it is assumed to have a finite number of trail- 
ing edge flaps, NF in number, which provide for the control torques. 

The flaps are assumed to function as rigid bodies which are good approx- 
imations for most control flaps . The wing is assumed to have a large 
aspect ratio, at least, for structural purposes. This permits the 
wing to be modeled as a beam-rod. The actual root is replaced by an 
effective root normal to the elastic axis of the wing. Rotational 
inertia, shear deformations and sectional bending are neglected. 

Let h(x, t), e(x, t), P^(t) i = 1, 2, 3, ••• WF, be the vertical 
displacement of the wing along the elastic axis, the torsional dis- 
placement of the wing along the elastic eixis and the angular displace- 
ments of the flaps respectively. The kinetic energy of the wing/flap 
system may be written as 
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1 

2 


J (m (x ) ( 
0 


c)t 


(x,t))'+ I 




(x. 




h (x ,t ) 


W 

E I i 

i=i 


't) )"+;:>(!, -I- ]);'■ 



f'(x^t) 


+ 



• 1 ) 


where is a.ssuin,:^ i to be equal to xo''o i. 'egions where the 

fl&p do0s not xxi-h-* 'i XG po i-GTix, I ctl v^i\- •;.• t r- .^y s tjerT'. may bo written 

as j, 


1 ^ 't^h(x-t) G ; 2 

- 2 I (EI(x)( — ^ i (tJ(x)i, ) )dx 

0 ' "" 




1 






A prime (') will be used to represent. xni :■. ^ to represent 

c) 

^ * 


The aerodynamic fer:'--:; and control L 'na,- v- riewed as non- 

conservative force;. - 

Using the .baprarrpan +xc.'hriiqir^ , c- Iciv- r^ .'^'' ’iovnng equations 
of motion and bountlary ^ ‘ond i t. i ^ns , 
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NF 

m(x)K(x,t)- S (x)e(x,t) - ( r Sg 0 (t))+(El(x)h"(x,t))" = L.(x) (2.5a) 

i=l " 

NF 

-S (x)K(x,t)+I (x)0(x,t)+( E (I. + b(c.- a)Sp, )P. (t))-(GJ(x)e ' (x,t)) ' 

i=l ^ ^ 

= Mg(x) (2. 3b) 

Kp 3.(t) 

-Sg fi(x,t)+(lg + b(c - a)Sg )0(x.,t)+ I 0 (t)+ — ^ 
i i ^ ^i" *i-i 


“ (x)"!' Tcg ^ i - l,2j3^. 

Pi Pi 


•NF 

(2.5c) 


where Lj^(x), M^(x) and Mg (x), i = 1, 2, 3, ... NF, are the aero- 
dynamic forces and moments and TCg ,i=l, 2, 3, ...NF, the control 
torques on the flaps . 

The boundary conditions at x = 0 are given by. 


(El(x)h"(Xjt)) Kj^h(x t) = 0 

(2.4a) 

El(x)h"(x,t)+ Kj^h'(x,t) = 0 

(2.4b) 

GJ(x)e’(x,t)+ K^e(x,t) = 0 

(2.4c) 


The boundary conditions at the wing tip x = iw are 

El(x)h"(x,t) = 0 (2.5a) 

(El(x)h"(x,t)) ' = 0 (2.5b) 
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GJ (x)e ' (x ,t) = 0 


(2.5c) 


If El(x) and GJ(x), the flexural and torsionsuL rigidities are 
zero at the wing tip, the boundary conditions are, 

El' (x)h"(x,t) = 0 (2.6a) 

2EI’ (x)h" ' (x,t) + EI"(x)h"(x,t) = 0 (2.6b) 

- GJ ' (x)e ' (x,t) = 0 (2.6c) 

In the above equations the inertia and aerodynamic loads are also 
assumed to be zero at the wing tip . 

Since the properties of a practical wing are not uniform it is 
difficult to solve the above boundary veQ.ue problem even in the 
absence of external loads . Hence these differential equations are 
reduced to integral equations and then discretized. 

The influence functions C„„(x, §), the static deflections at 
station x due to a unit load at station § along the elastic axis , and 
Cgg(x, I ), the state torsional deflection at station x due to a tmit 
moment at station I along the elastic axis, for a wing with a flexible 
root, are given in Appendix A. 

The inertial properties are assumed to be lumped at NB collocation 
points, x^, along the elastic axis. Then following standard techniques 
it is possible to rewrite the equations of motion, in the absence of 
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external loads, in the form. 


[l]fqo 3 - -[Kq] [MQ]{q^) ^2.7) 

where 

[l] is the unit matrix, 

is the vector of nodal displacements, torsional 
rotations and flap rotations , 

[Kq] ^ is the Influence coefficient matrix, and [Mq] is 
the inertia matrix. 

Defining [U ]'^ = [Tc„ , Tcq , . . . Tc. ] 

1 ^2 

m KF X nB • NF X NB • NF X NT 
and [Bq] = [ 0 : 0 : I ] 

we may rewrite the entire eqmtions of motion as 

[K]_]q]_ + [NLj_] = F^Ct) + BU^ (2.10) 

^ere F^(t) is the aerodynamic load vector, is the vector of con- 
trol torques, B the control torque distribution matrix, and q^ the 
vector of modal displacements . It is also assmed it is possible to 
measure modal displacements and rates and these measurements can be 
written as. 


( 2 . 8 ) 

(2.9) 


= %qi + H^q^^ 
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If a linear aerodynamic theory is assumed, the Laplace transform 
of this load vector, Fq(p), vdiere p is the Laplace transform variable, 
is specified as , 


Fo(p) = I p„A^-b-[Q(s,M)]{q^(p)} 

where Q,(s,M) = {q. .(s,m)} is an aerodynamic load matrix, defined by 
the relations , 


qij(s,M) 


1 J.J. 

" \ 2 Pi v" 


where is the area of the lifting surface for a three dimensional 

lifting surface and equal to the chord for two dimensional aerofoils^ 

z^(x,y) is the modal deflection s\irface in the i^^ mode and Ap'^ (x,y^ 

th 

s,M)^ is the Laplace transform of the pressure difference in the j 
mode . The modal deflection iSurfaces may be obtained by solving 
the free vibration problem given by equation 2*7 • 

A general theory for approximating the aerodynamic loads q..(s,M) 
in the frequency domain, will be presented in the next chapter. This 
theory will permit the entire a^roelastic system to he replaced by a 
finite state modal . 
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i'lNr'Ii; STATK MODKI.ING Oi' AKKODYNAMTC l.OADS 


In thi;.; fieri ion, a general method of approximating the aero- 
dynamic loadc in the frequency <iomain is presented for wings in two 
and thx’ee dimensional llov;. 

Twi) 1) i liu-i 1 j ( ma 1 AirftM Is 

'L’he i^ade oxlri-ant nx't.hud. Is applied t,o cbiain the solution 
to the- prossure hi la-cributlon nn wings perlonrsing arbitrary oscillations 
in two dime ns -i.o rial ble llov;. The int-thod is then generalized 

and applied to tv/o (limens i anal convoressibie flov; pr‘Ob"^ems 

^ i b 1 i 1 ow , -• Aerodynnin i c 1 oiuis (>ii aT*t>i t rari 1 y 

oscillating rigid v/iiiy,s in two-d i mens i ona 1 inccnipress i h 1 e ran be syn 

thesized from Wagner- hj solution for the li ft on awing due t.o a sudden step- 
wise change in the downvvarh as described in Kef. 1. R. T* Jones [10] 
has obtained an approximate Laplace transfoimi of Wagner Indicial 
function. Jt can be shown that eouations lot the pressure distribution 
and aerodynamic loads may be obt.aixned for’ convfn:‘glng or diverging 
oscillations by rc: placing 0 (k ) by b(s) s<fAs) where p(s) is the 
Laplace transform of Wagner's Indicial function. Thus equations of 
motion of a wing with a trailing edge flap (id] or a wing with an 
a. lei: on anl a. 1 1 ay) [. y. . W may br; obtained by using the correspond- 
ing equations for- oser 11 1 at:ory motion. T^or converging oscillations 
the wake is assumed to be finite and dependent on initial conditions . 

In this sec lion Theodorsen's c*irculation lag function is 
anaJytically continued f'or conv' 'rg;i n,g i--n nt' cci cirfei] wifdi no 
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oscillating components in two dimensional incompressible flow, and 
the physical significance of the results explained. 3y representing 
Theodorsen’s function as a series of Rummer f-unctions, an asymptotic 
expression for Wagner's indicial function is obtained. The first two 
terms of this approximation are identical to the approximation for 
Wagner's function obtaired by Garrick. By applying the Fade 
approximant theory it is shown that Theodorsen's function may be 
represented by a sequence of rational functions which converge 
uniformly. This method of approximating Theodorsen's function is 
generalized so one can construct rational function approximations not 
only to Theodorsen's function, but to gust response problems also. 
These rational function approximations may be easily inverted, using 
Laplace inversion, to obtain good approximations to various indicial 
functions . 

The lift and moment about the elastic axis on a rigid thin air- 
foil, performing vertical translational or torsional oscillations can 
be easily obtained by integrating the pressure and its moments (given 
in Appendix B) . 

These expressions are, 


r 1 


Lb 

^ -l2 

M 

O' 

= 27Tqb 


s + 2si/)(s) 


s - as^+ 2ip(s ) ( 1 + ( i - a)s ) 


as +2(a + -)0(s) -s |s (a - g)-(^+ a‘~)s"^+ (a + ^)2-ip(s) (1 + - a)j 


X 


h/b 


o 
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These equations may be written also as , 


r 1 

(r 

Lb 

. ,2 1 

M 

= 2mqb i 

a 

/ 

* 

u 


2 

s 

2 

as 


2 

s - as 

sCa-l/2)-(l/8 + a^)s^ 



2 


" h 

+ 


0(s) [s, (l+( - - a)s) ] ) 

b 


2a + 1 

1 

or 


- 

) 

- 


Thus the circulatory part of the aerodynamic matrix can be decomposed 
into the product of two vectors • As a consequence we have ^ 


Lb 

p 

s^ s - as^ 

or 

= 2TTrib 

as^ s(a - |)-(^ + a^)s^ 


h 

b 

a 


+ 2Tiq,b^ 


2 

2a + 1 


y(s) 


y(s) = ^(s)[s,(l +( I - a)s)] 


h 

b 

a 


Thus it is sufficient to construct a finite state model of ^(s) (shown 
in Figures 2 and 5 ) in order to obtain finite state models of the 
generalized aerodynamic forces • An application of this approach to 
flutter suppression studies is presented in Ref. 30. 
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In general; if the influence of compressibility and finite span 
are included; it is impossible to decompose the matrix in the above 
manner* In two dimensional incompie ssibls riow, the circulatory aero- 
dynamic forces are dependent only on the iowavash velocity at one 
point (the chord point for a rigid airfnil ). In compressible 
flow or for three dimensional lifting surface.''; the aerodynamic forces 
are dependent on the downwash at a finite uumber of collocation 
stations. Further the contribution of the iO'wnwash at each of the 
collocation point; to different generalize! forces, is different. Thu 
one cannot assume apriori that a decornpositi on of the form discussed 
above is possible. Thus, the above epp; 't is not useful in com- 
pressible flow or for three dimensional lifting sui'faces . 

Garrick has shown that Wagner indicial function and 
Theodorsen’s function are related by the Fou>-i.er transform relation- 
ships . Thus ; 


%(t) 


1 + i r 

2tt - ik' 


ilvt 


for t > 0 


0 for t < 0. 


This may be written in the form 

oo 

Wg(t) K^(t + T)- 1 - f (if 


where ( ^ ) 


C(k)- 1 iKT 
ik ^ 


_00 



This equation suggests that numerical quadra tiire may be used to 
compute K(t) for T < t < 2T for any given T. An equivalent possibility 
is the use of discrete Fourier transforms like the fast Fourier trans- 
forms [ 51 ]' In order to show the disadvantages of this technique a 
brief discussion of the method is presented. 

To observe the effect of sampling both in the frequency and time 
domain, at finite but equal intervals, the above equation is evaluated 
at the points , 


t = J-At, j = 0, + 1, + 2, ..., + (N-1) and 
J 


n 


n=0, +L,+2, (N-1) 


where N = TF^ with F = l/At and T = l/Af- At. is chosen small enough 

j 

for the highest frequency present to be sampled at least twice dviring 
each cycle. With the assumption that the indicia! function and the 
spectrum C^(f) are periodic with period T we have. 


N-1 


W (d-At) = ^ 2 C (n.At)e’''^'^ 

C,(n-Af) = |''E\(d.At)e-2^ 
rt=0 ® 

Let A = C~(n,Af), X. = TW (j,At) and w = 

J 8 


We then have the discrete transform pair, 

N-1 

X. = s A^w ‘’fx.3 = [n]{A^} . 
rt=0 ^ 


N-1 ""J 
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From the orthogonality properties and the properties of the exponential 
function it is possible to factor [fi] in such a manner that the number 
of computations involved for inversions are reduced from to ^ log^K 
(provided N = 2’^). This decomposition technique is the Fast Fourier trans- 
form. Thus it is possible to obtain numerical values for W (t . ) for 

6 o 

various values of t. by using different sampling periods T to any desired 
degree of accuracy. Also errors due to the assumption of periodicity 
can be corrected for^using Window Techniques. However this does not lead 
to a finite state model • From the numerical data generated for the 
indicial functions;, K^(t), it is essential to fit the data with a finite 
number of exponentials (2 or 3 )> This is usually a non-linear technique 
[52]. The Laplace transform of such an approximation would then lead to 
a finite state model in the frequency domain- This whole procedure is 
xindoubtedly extremely tedious and involves several computational steps . 
Conceptually the procedure involves approximating a distributed parameter 
system by a discrete time system which is then approximated again by a 
continuous system in the time domain. Thus a direct procedure of finite 
state modelling from known values of the response at f = f^, with the 
correct initial and steady state behavior is desirable . In order to develop 
such a technique, Theodorsen's function is analyzed in detail - 

Analytic continuation of Theodorsen's function: The function i(j(s) 

will now be evaluated in the left half of the 's' plane along the real 
axis and the physical significance of the results will be explained. In 
this case, ip(s) may be written as 
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tj)(s) = 1 - 


1 


^ r 

p^-s cosh V 

1 + Lt cosh V dv 



f -s cosh V , 

J e dv 

0 

The integrands of the two integrals are monotonlcally increasing in the 
domain of v. Further they are both divergent. Hence the limit may be 
evaluated by using L 'Hospitals rule, which immediately leads to 

= 1 for Re (s) < 0, Im(s) = 0 . 


This result is of fundamental importance . It was assumed that at time 
t = 0, the wing was stationary and unperturbed. Clearly, converging 
motion with no oscillatory components from an unperturbed state is not 
possible. This implies that on the negative real axis of the 's’ 
plane the 'motion' remains unperturbed. As 0(s) ^ 1, is the solution 
for the case of steady translation of the wing, the quasi-steady solution 
for the pressure distribution is the exact solution on the negative 
real axis of the 's' plane. 


Asymptotic expansions of Wagner's Indlcial function: In terms of 

Kruramer functions we have 


t(s) = 

1, 2s) 

sU(li, 5, 2s) 


Re (s) > 0 


The Kummer functions U(a, n, z) are a class of confluent hypergeo- 
metric functions that are multiple valued with the principal branch 
defined for -v < arg(z) < tt. The order of the singularities for large 
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and small values of z are determined by a and n respectively. 

As a consequence of the properties of Kummer fiinctions ^(s) may be 
represented as an asymptotic series for large s as^ 

t(s) = l -^(4su(l,l.^s)+ E B (C.s)u(l,l,C.s)+ 0( i) 

i^l 1 X IS 

N = 0, 1, 2, ~ . 

where and are constajits • Apart from being asymptotic for large 
s the expression also satisfies the condition, 

Lt t(s) = 1 . 
s-^0 


Further, the inverse Laplace transform of the asymptotic series can be 
easily found. Hence it is also possible to find an asymptotic expression 
for Wagner's indicial function. In particular, for W = 0, 



where P(t) = 1 + 0.875t + 1.281^55t^ + 1.8l^28>^ + 4.09134 t^ . 


The first two terms are identical to the indicial fxinction given by 
Garrick. In general the correction to Garrick's approximation may be 
written as , 



(Polynomial of degree N in t) 
(Polynomial of degree N+4 in t ] 


The rational function in the brackets is usually referred to as a 
[N + 4 , N] Pade approximant. Baker [ 33 ] has discussed the methods of 
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construction and the properties of Fade approximants . 


Pad4 approximants of Theodorsen's function: For small values of 

W. P. Jones [J4] gives the following series fort|r(s): 

^(s) = 1 + ys + y^s + s^( + 3^ ) 

where y = log (s/2) + Y- The correct asymptotic series for small value 
of s is given by. 


t(s) = 1 + ys + y^s^ +s^[^-^+y5+y ] + s^(y^- y^+ 


.2 


y - 


3 4, 5 5 19 2^ 11 




+ S^(y6. 2y5p ^y^ ^y5p ^y2- gyp O(s^), 


Re (s ) >0 . 


Obviously it is quite difficult to construct a Fade approximant from 
this series. However for large values of s ^ we have the asymptotic 
expansion^ 


t(s) = |(1 + ^ - J:^ + T _ ^ _ __6^ ^ 

8s 64s^ 128s 512 s^ 1024s 4194504s' 

Re (s) >0 . 


The fact that the lag ftmction has a value of 0-5 for s = “ and 1 for 
s = 0, suggests that it may be approximated by an [N,N] Fade 
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approximant quite accurately- To demonstrate this, the first four 

Fade approximants are constructed from the asymtotic series for large 

s. They are asymptotic to the lag function up to the first, third, 

fifth and seventh orders in — respectively. They are given by 

s 


s + 0.5 
2s + 0.5 


+ l.^s + 0-37^ 

2 ^ + 2.5s + 0.575 


s^ + 5-5s^ + 2 .7125s + 0.^4-6875 
2s^ + 6.53^ + 4.25s + 0.46875 


and 


s^. + 4.64696s^ + 9-55571s^ + 5-5l755s + 0.49554 „ / > ^ 

— r 5 ’ P , Re(sj > 0 , 

2s + 8.79592s'^ + 16.71894s + 7.67296s + 0.49554 


These expressions are found to converge rapidly over the entire right 
half of the plane. However the above approximations are not very 
good for small values of s as they are obtained from the asymptotic 
series for large s . 

This method of approximating Theodoresen ' s function and evaluating 
the Wagner indicial function approximately is similar in principle to 
the method suggested by Luke [55] for evaluating the Randall function 
in the theory for the unsteady aerodynamics of oscillating cylindrical 
shells in supersonic flow. 

The method of construction of the Fade approximants may be 
generalized. Since the lag function is known exactly along the 
imaginary axis s = ik^ the Fade approximant may be constructed by a 
suitable least squares technique. If C(k) - N(ik)(D(ik)) where 
N(s) and D(s) are polynomials in s^ the coefficients of tne polynomials 
may be determined by minimizing 
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s (D(ik.)C(k,)- N(ik,))^ 
j a 0 


where kj,J = l, 2 ^ 5 ^ ...^m are the vaJ-ues of k at n different points- 
This technique needs only the values of C(k) at a finite number of 
points in the frequency domain^ and hence may be generalized to other 
cases ^ where the aerodynamic loads can be calculated numerically only. 
This problem may then be reduced to the solution of an over determined 
linear system for the coefficients of the polynomials N(s) and D(s). 

This procedure results in 

^ sV 0.7610565^+ 0.1020^8s^+ 0 . 002^^0673 + 9. ^^732 X 10'^ 

2s + 1.0659963^ 0-115928 s^+ 0.00261680s + 9.55752 X 10"^ 

Re(s) > 0. 

Clearly Wagner's indicial fimction may be approximated using the 
formula^ 

K(t) = = 1 - E A.e 

" i=l ^ 

where are the roots of the denominator polynomial where 


i 

A. 

1 

^i 

1 

0.011285763 

0.0044482234 

2 

0.043280564 

0.027697193 

5 

0.21639860 

0.096054968 

4 

0.22903508 

0 . 40379780 
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ESy taking the Laplace transform of the approximation for Wagner's 
function obtained by R. T. Jones [10], we have 


ip(s) 


0.^6l^s + 0.0273 

2s^+ 0.6910s + 0.0273 


Re(s) > 0 


The [4^4] Fade approximant obtained by least squares and the corres- 
ponding indicial function are compared to the exact values and the 
results of R. T. Jones (Figures 4-8). 

This approximation is very good both in the frequency and the time 
domain. The simplicity of this method compared to the one given by 
R. T. Jones is quite obvious. The main advantage of this method is that 
it can be generalized to three dimensional lifting surfaces- 

It is also interesting to note that the poles of the six approxi- 
mations for the lag f\inction are negative real values - This indicates 
that the transient pressure response to a downwash input is 
asymptotically stable . Further since the negative real axis is a 
branch cut of the exact expression for the lag function it represents 
a continuous distribution of poles- 

Application to gust response problems: In applying the above 

theory to gust response problems, a certain amount of caution is essential. 
It must be emphasized that only the Laplace transform of the temporal 
portion of the circulatory pressure distribution can be approximated by 
a sequence of [n,N^ Fade approximants , This does not hold for sinu- 
sodially convected gusts. 

Sears [9] considers a gust velocity of the form^ 
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w - w 

'^G " '^0® 

€uui shows that lift to be equal to 

L(k) = k-ir q-b- e^"* S(k) 
where 9(k) = (jQ(k)- ij^(k))c(k)+ ij^(k). 

Kiissner's indicial function is related to S(k)e in the same manner as 
Wagner's indicial fimction is related to C(k). However to obtain an 
exponential approximation of Kussner's function one shovild construct 
a [n,N- 1] sequence of Fade approximants of S(k)e”^^. In fact, using the 
least squares smoothing technique developed above and taking the inverse 
Laplace transform we obtain the following expression for Kussner's 
indicial function K^Ct), 

^ -0. T 

K2(t) = 1 - A.e ^ 


i 

A. 

1 

Pi 

1 

0.012994467 

0.0049896174 

2 

0.062519920 

0.034931400 

3 

0.40920539 

0.13796713 

4 

0.51548022 

1.1645811 
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In order to show the relationship between the temporal and spatial 
portions of the gust response, it is essential to consider a gust of 
form, 

W = W e^"^ " 

G 0 

Kemp [ 36 ] shows the lift to be of the form 

w 

L(k,X) = 4 ttxi b e^“^^ K(k,\) 

where K(k,X) = ((Jq(X)- iJ^(X))c(k)+ ^ Jj^(X)) 

and X = . 

The first term in the function K(k,X) is the circulatory part while 
the second is the virtual inertia effect. Clearly it is now sufficient 
to approximate C(k) by an [N,N] Fade approximant to obtain the indicial 
function corresponding to a downwash of the form 

H(t) 

where H(t) is a step function. 

A sinusoidally convected chordwi.se gust past an airfoil at a 
steady angle of attack a, of the form^ U = known to 

u U 

produce a lift equal to [ 57 ], 

L(k) = liTTq abe^ T(k) 
where T(k) = C (k ) ( (k ) -ij^ (k ) )+ J^Ck)! i2J^(k). 
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From this result we may obtain the indicial function for a sharp edged 
chordinse gust striking the leading edge of the aerofoil at t = 0 


U X > Vt - b 

Uq X < Vt - b 

The indicial lift is given by 

L(t) = 4TTq ^ «bK^(t) 

where K^(t) may be found in the same manner as Kussner's function. In 
fact K^(t) is approximately given by 

4 -0 . T 

K,(t) = 2 - E A. e ^ 

5 ,• 1 1 



i 

^i 

Pi 

1 

0.026965564 

0.0023690134 

2 

0.044268810 

0 .o 4 o 647302 

5 

0 . 45282297 

0.14614369 

4 

1.4759427 

1.7305802 


Separating the spatial and temporal portions of the gusts into 


U 

g 


iwt - ivx/v 
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we may show the lift to be 


L = Uim obe^^^ V(fc,X) 

where V(k,\) = C(k)(jQ(X)- iJj^(X))+ Jq(X)+ iJ3^(X)+ 


~X 


Sharp edge gusts moving with respect to inertial observers are 
considered} i.e-^ gusts of the form 


and 


( 0 X > Ut - b 

/ Wq X < Ut - b 

0 X > Ut - b 

Uq X < Ut - b 



where U is different from the free stream velocity V- is asstnned to 
be very small compared to U so only first order effects need be con- 
sidered. The corresponding indicial functions are related to K(k,ak) 

e and V(k^ori)e ^ where a = V/U, in the same manner as Kiissner’s 

-ik 

function is related to S(k)e and hence may be written as_, 

^ T 

= 1 - E a / e ^ 

J=1 

“ c 

These indicial functions are shown graphically in Figure 9 and Figure 

10 . 

Before going into the natvire of the solution of the pressure distri- 
bution in two dimensional compressible flow it may be instructive to 


32 



review the effects of thickness and viscosity on 0(s). Hewson-Browne 
[38] has investigated the oscillation of a thick airfoil in an incom- 
pressible flow and obtain a general expression for ^(s) (which he 
evaluates only for simple harmonic motion I) including thickness effects. 
The generalization is different in the expressions for lift and moment. 
This would destroy the decomposition property of the circulatory portion 
of the aerodynamic matrix discussed earlier. The singularity in the 
derivative of ^(s) continues to exist. Woods [39] has presented a new 
approach for the calculation of the unsteady two dimensional flow about 
aerofoils performing arbitrary motion in an incompressible fluid. While 
the flow is assmed inviscid^ the potential flow boundary conditions 
are modified semi -empirically to make some allowance for viscous 
effects. The method is applicable to thick airfoils, the only limita- 
tion being that the velocities and displacments of the 'unsteady pertxirba- 
tion about the mean motion be small. Results are obtained for the lift 
and moment by an application of Blasius theorem. However these results^ 
while being very impressive mathematically seem to be practically useless 
as most of the integrals cannot be evaluated for general motion. Pro- 
vided the Reynolds number is sufficiently large for boundary layer 
theory to be applicable, viscosity has three main effects on the 
theoretical potential flow. These are, l) the Kutta condition is 
modified in that the position of the rear stagnation point is indepen- 
dent of incidence 2) the velocity distribution of the mean steady flow 
particularly near the rear trailing edge is modified 3) viscosity 
contributes to the damping that is independent of the velocity but does 
not effect the inertial loading. Woods accounts for 1 ) and 2 ) but not 
3)* The results indicate that the reduction in the wake velocity due 
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to the thickness and viscous effects^ changes the flat plate deriva- 
tives by an ajnoimt which is quite large for k large but neglegible 
for small values of k. Chen and Wirtz [^0] have included second order 
wake effects^ but it appears that the singularities they obtain are 
no different from those of the first order solution. 

It seems probable that the nonlinear and viscous effects in the 
wake may reduce the wake velocity and distort the wake which may be 
significant especially for low values of k and large values of [s]. 

So far there is no theoritical justification of this* The above consid- 
erations are important if experimental verifications of Wagner’s solu- 
tion are sought. 

Subsonic flow . - In Appendix C, the behavior of the pressure 
differential across the airfoil in subsonic flow was shown to be 
essentially like incompressible flow for low values of reduced frequency 
and mach number. However for large values of reduced freq^iency^, unlike 
incompressible flow^ the force required to generate impiilsive motion may 
be expected to remain bounded for all time and finally reach a steady 
state value. In particular^ the work required to generate impulsive 
motion is not entirely recoverable^ so that the process is irreversible^ 
and some part of the starting force must be regarded as ’’damping’’. It 
is well known ^ that for large values of reduced frequency the pressure 
is given by Piston theory [1] for all mach numbers, M, greater than 



idiere z is the displacement of the airfoil. 

Another curious fact, that has not received the attention of any 
previous investigators in unsteady aerodynamics, is that the steady state 
aerodynamic load on an airfoil in steady subsonic flow or supersonic 
flow is always finite as long as the mode shapes are bounded. This 
indicates that the aerodynamic system by itself is a completely stable 
system, in subsonic or supersonic flow, in spite of being a circulatory 
system. This physical fact fully explains the reason for all the poles, 
of each one of the Fade approximants of t(s), being in the left half of 
the 's' plane. The stability problem is especially significant in 
transonic shock free flows and hence is discussed in greater detail in 
the next section. 

In Appendix B, it was shown that the integral equation for the 
pressure distribution can be reduced to a Fredholm integral equation of 
the second kind with a completely continuous kernel. The significance 
of complete continuity will be discussed in 5.2. In Appendix D, it was 
shown that for an eqxiation of this type, the solution has the form of 
a [N, Iff J ] Fade approximant. Clearly, comparing this result, with the 
Fiston theory result for large s , we conclvide that the loads must be of 
the form. 


q. . (s,M) = It 
N-»“ 


^Ntl, .N, 

V +v+--*^i 


„ ,N, - N-1 

Ms + bgS 


.. b. 


m-1 


■vdiere the coefficients a^ and b^ are functions of Mach number alone 
Clearly as s ■-><», q^ . = 8.^/ii and a^ can be obtained from piston 


theory. In the limit as s 0, q^^. (0,M) = (M) and 

equals the steady state value- As the mach nvimber approaches zero. 
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one of the poles of the denominator moves to negative infinity in the 
's’ plane. This clearly indicates that in incompressible flow an impul- 
sive force is needed to generate impulsive motion. Hence for large values 
of s, q_ behaves like an inertial force in incompressible flow. In 
compressible flow, unlike incompressible flow, the coefficients b^ are 
different for different types of loads and mode shapes . 

This difference in behavior may be easily explained physically. 

In inccanpressible flow there is a lag in the development of the circulation 
around the aerofoil, while the pressure variations are propogated with an 
infinite speed. Incompressible flow in addition to this lag the 
pressure variations caused by the motion of the aerofoil are propogated 
by a finite speed. This introduces an additional lag in the loads. 

On the other hand as the speed of the airfoil increases, it is clear 
that the pressure variations move upstream more and more slowly in 
relation to the aerofoil . If the airfoil moves faster than the speed 
of sound, it cannot make its presence known to the fluid ahead of it and 
the mechanism for bending the streamlines ahead of the airfoil no longer 
exists. Eventually for very high speeds, the influence of the pressure 
variations is only local. At this stage the airfoil behaves like a one 
dimensional piston and piston theory leads to the exact pressure on the 
aerofoil for all values of s . 

Thus 

aQS 

lij(s,M)= for M “ 

It was also shown in Appendix D that the sequence of rational 
functions for increasing values of N converge and hence can be 
truncated and assumed to be of the form 
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(3-1) 


(s ,M) = 


mi, N 


0 


+ V + ‘iij (">Vi 


„ N ^ , N-1 , , 

Ms + bgS . . . + 


for some fixed N« 

Further, in view of the aerodynamic system being completely stable, 
we may conclude that all the poles of q. ,( s,M) lie in the left half 

j- j 

of the 's' plane. This is an important and useful property, as it can 
check on the validity of the rational function approximants after they 
have been constructed. 

The exact solution for subsonic flows , can be obtained by Mathieu 
functions [1]. However the actiial computation of the solutions is ex- 
tremely tedious. Firrthermore , the solutions are not closed form 
solutions and hence are only approximate. Thiis it is difficult to 
obtain series expansions in frequency of the aerodynamic loads, in order 

to analytically construct the Fade approximants. 

Hence a numerical technique used often in the area of system iden- 
tification and for digital filter synthesis [4l] is proposed and 
applied to various cases • The problem of computing the approximants 
is reduced to a least sqviares minimization problem. This is then 
transformed to a overdetermined linear system for the coefficients of 
the Fade approximants 8ind the solution obtained by QR decomposition. 

Now q. .(s,M) may be easily calculated by the nmerical method out- 
lined in Appendix C for a given mach mmiber, and for various real values 
of reduced frequency. 

k = k. = -is m = 1, 2, • • ' > ’ 

mm ^ 
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Also Bq can be obtained from piston theory and q5‘. (m) from steady state 
calculations • Thus the above equation may be written as , 




-1. -2 , -N 

s + s a_ + . . . s a 
m m2 m N 


- q 


10 


(s ,M)(b^s“‘' 
m^ ^ ^ 2 m 


V"’') 

K m ' 




'0 


m = 1, 2 , 5, ... m, 


0 


(3.2) 


Hence these equations are a system of complex equations in the 2N real 
unknowns a^.^ a^, ... a^^ and '^ 2 .’^^ ^rj+l' ”"‘0^ chosen larger 
than N, this system of equations may be solved by the linear least 
squares technique . 

Numerical experience indicates that a large number of values must 
be chosen for k^^, in the low and high freqtiency regions . In the low 
frequency regions the coefficient of could be a small complex 

quantity while at high frequencies the right side of equation (5.2) 
would be a small complex quantity. Hence a large number of values 
must be chosen for k^, in the low and high frequency regions, in order 
that the approximation be sufficiently accurate in these frequency 
regions . 

These approximants for nondimensional generalized forces in (l) 
plunging, (2) pitdiing and (5) flapping mode and tabulated in Tables 
1 "* 3 • The flap hrnge xs assumed to be at the rear quarter chord point. 

In most of these cases N was chosen equal to 2 . In these cases it was 
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found that more accurate results could be obtained with N = 4. However 
these results are not tabulated. 

It was assumed that it is possible to calciiLate the aerodynamic 
loads q.. ^(s ,M) acciirately for all values of k and M in the subsonic 
regime . This assumption is generally true in two dimensional subsonic 
flow, provided a fairly large number of collocation and integration 
points are used in the method described in Appendix C- The method is 
very accurate even in the presence of control surfaces, provided 
collocation points are not chosen too close to the hinge line . 

The influence of compressibility which is well known in steady 
flow is different for the steady and unsteady components of the loads . 
For low mach numbers (O.l < M < O.J) it was found that higher order 
Fade approximants were required to model the pitching moments , where 
the poles of the approximants have a tendency to move towards the 
imaginary axis of the s plane • This could be due to the higher order 
of the approximants without any physical significance. 

For high subsonic mach numbers ( 0.7 < M < O. 9 ) the loads were less 
stable even with lower order Fade approximants • On the other hand it 
is well known [42, 45 and Section 9*2, 1] that oscillations of a two 
dimensional aerofoil about a pitching axis may become unstable for low 
values of reduced frequency and high mach numbers, if the axis lies in 
a certain region, ahead of the quarter chord axis and if the moment 
of inertia about this axis is sufficiently large. This has been 
attributed to the fact that the logarithmic term in frequency in the 
kernel fimction is stronger at hi^ mach numbers for two-dimensional 
aerofoils than for three dimensional lifting surfaces. 

Hy taking the Lapalce inverse transform of q. .(s,M)/s, it is 
possible to obtain the indicial response in the i^h mode due to a time- 
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4*V^ 

wise step change in displacement of the j mode at time t = 0. These 
results are shown in Figures 11 to 19, and have been obtained from 
[1+, 5] Fade approximants in all cases except in the case of at M = 0.7* 
From these figures it is quite clear the loads corresponding to the same 
pressure mode have the same general behavior in the time domain. Also 
the loads corresponding to the flap mode take much longer to reach the 
steady state value, than the loads due to plunging and pitching. 

The influence of compressibility on Wagner's indicial function may 
be studied. In this case, the lift q-j^-]^ , due to a rigid pliinging mode is 
calculated at several Mach numbers and the corresponding indicial 
function is obtained by taking the Laplace inverse of • This 

result is shown in Figure 20. These curves are compared with that of 
Mazelsky and Drischler [15] for Mach number :::: 0.5, in which case the 
error was found to be the highest. 

Supersonic and transonic flow . — In supersonic flow the 
expressions for the lift, moment and partial moment for an aerofoil 

with a flap are given in [4k] and for an aerofoil with both trailing 
and leading edge flaps in [45]. The Fade approximants for the exact 
expressions for the lift and moment of a rigid aerofoil in plunging and 
pitching may be obtained. The lift coefficient for an airfoil in pure 
plunging motion is , 

L, 2 

h.,, s + as + b \ h 
_ = 4( ) - or 

Ms + cs + pb 
4 / s + Pm - s h 

MS + . ^3 N 
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where 


+ 5 + JM-'O-M)- A = M^(7M^-5) 

2 2 

b = -^ (4lV(P -M)+ 3P^), c = ((3M^+ 1)(P -M)+ IJM) 

The moment coefficient in piire plvinging motion is , 

^ = 4( — ^ , where ^ 

qb Ms + c^s + b^ 5M(11M - 5) 


5P 


^ M(llM^-5) 


and 


lOMp 
IIM^- 3 


The lift coefficient in pitching is given by. 


L 


a 

qb 




a^s + b 
Ms + pb 


2 

2 


)a , 


, M. 3M - 2 X , 

' 2 ■ p' ;at ’’ 2 

3M + 1 


3MP 

3M^+ 1 


The moinent coefficient in pitching is given by, 


2 

M s + a^s 

-2L 4f ^ 

^2 ^ 5Ms + b. 

qb 5 


2 

)u , a ^ - p ^ ^ ^ 


3(2 - m^)mP^+ 3P^ 


2M 


All four approximants have the correct asymptotic behavior in the 
limit s ^ 0 and the piston theory limit was enforced for s = Also 
the poles of all four of the approximants satisfy the stability cri- 
terion discussed in the last section. The approximants are good 


41 


approximations for the loads for M >2.0. Thoiigh not extremely accur- 
ate these approximations are certainly better than piston theory or 
quasi-steady theory . 

It was observed that for higher order Fade approximants , the poles 
behaved in an erratic manner especially for low supersonic mach numbers, 
moving often into the right half of the ^s' plane. In transonic flow, 
it is well known from the exact solution for the pressure [46], that 
the linear aerodynamic theory leads to infinite aerodynamic loads in 
the steady state as the loads have a pole of the form /s • This was 
probably the reason for one of the poles of the approximants moving 
towards the origin for low mach numbers in the supersonic regime and 
high mach numbers ( > 0.8) in the subsonic regime. For the reason 
indicated above it seems appropriate that in transonic flow, the 
aerodynamic loads may be written as a [N,W2] Fade approximant in -/s 
with a pole at ^s = 0 . However this result has no practical value , as the 
linear theory is known to be invalid in steady transonic flow. 

The physical interpretation of this result is that two dimensional 
shock-free transonic flows are unstable on the basis of linear theory. 
One of the arguments against the physical existence of shock-free 
transonic flows was based on the supposed instability of these flows 
with respect to upstream moving disturbance waves [47]. The argument 
was that these waves when superposed on a steady shock-free basic flow 
could move upstream as long as the local steady flow speed was subsonic, 
but as they entered a region of supersonic flow they would come to a 
standstill and coalesce to form a steady shock wave. 
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However Pearcey [ 48 ] showed using laborious experimental tech- 
niques^ that the disturbance waves moving upstream can penetrate the 
supersonic region, and that in fact, shock-free transonic flow 
could be experimentally realized. This is attributed to the non- 
linear turning effect of the waves moving upstream into the supersonic 
region, which was not theoretically accounted for earlier [49]- Then 
based on geometrical acoustics one may show that disturbance waves 
moving upstream can penetrate the supersonic region. The problem of 
shock-free transonic flow is too complex at present for a strict 
mathematical treatment . 

The numerical technique developed in the last section is also 
applicable to supersonic flow. The resialts are tabulated as in the 
subsonic case for N = 2 in Tables 6-8. It may be noted that the 
derivatives of the indicial response functions to time -wise step 
changes in the downwash in plunging and pitching [12] are not continuous, 
since the period of influence of an impulsive source is finite- The 
above technique approximates these derivatives by continuous exponential 
functions (Figure 21). Thus hi^er order Fade approximants are 
essential to approximate the derivatives of the above indicial 
functions . 

It was found that for Mach numbers less than 1-5 second order Fade 
approximants were insufficient for modelling the loads. The possible 
reason for this has already been discussed. On the other hand for 
Mach numbers greater than 2, a first order Fade approximant was sufficient 
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for the lift and moment in plunging and pitching in the frequency domain* 
This is to he expected as for a wide range of frequencies. In fact 
Piston theory may be considered as a zeroth order Fade approximant 
for high mach numbers. 

Three Dimensional Lifting Surfaces 

In this section a systematic method of obtaining Fade approximants 
of the aerodynamic pressure loads is presented. Though the method is 
actually quite difficult to apply in practice it will prove the 
existence of Fade approximants and the form of the solution. The method 
is also applicable for two dimensional airfoil s. 

First it must be noted that Piston theory is also valid for three 
dimensional lifting surfaces for high frequencies and that any solution 
for the pressure distribution must converge to this limit for high 
frequencies . 

In subsonic flow for two-dimensional airfoils, it is well known 
that the non-dimensional pressure Cp(^^Tl;,ik) may be written as ^ 

^ P 2^0 ^21 ^ ^ 5 - 5 ) 

For three-dimensional lifting surfaces it is given by 


c - ^ p"(ik)’^+ (ik)hog ik (ik)^ (3.4) 

^ rv=0 


where p^ and p^ are independent of f>-equency. Assuming the downwash 
to be Wqq + ik w^^ and using the kernel function expansion in Appendix 
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B, for low freqiiency (equation Bl8) , we have 


8tt(Wiq+ ikw^j^) 


= ]■]■[ E a (ik)^(ik)^log ik b (ik)'^]CL(C,T1,ik)dedi^, 

A ifcO " n P 


Hence we have the infinite system of equations , 

1 n' 0 ^ 

T' ^ 

0 = J'J'tVjo + V3°)‘>S 'll 

0 = + ^^50 + 

A 

The first equation in this infinite set, involves the calculation of the 
pressure distribution by application of steady state lifting surface 
theory. After having determined the steady state pressure the second 
term is also determined from steady state lifting svirface theory since 
p^ is known and p^ occurs in combination with the steady state kernel 
Eq. Thus the coefficients may be obtained numerically by a procedure 
like the kernel function technique. In this meuiner it is possible to 
obtain all the coefficients p^q and P^* Oiie special case of importance 
is when w^q =0. In this case p^ = 0 and p^ = 0. Hence in this case, 

^ = PjQ ik - PjQ k^ + O(k^) + 0(k^log ik) 


45 


Thus the existence of a solution of the form given by (5*2) is estab- 
lished . 

Next we note that 


s log{s^{s/sQ)) 


s 2 

i=l 



i 



+ s log Sq 


0 < s < 2 sq. 


From (5 '2) we have 


- 1 ) 


TT. n xn 


+ ik' 


- E (-1) 
n=l 


j+1 0 


^31 ( j 


This can be rearranged in the form, 

__ 2N-2 

C = s c^(ik)“ (3.5) 

n^O 

The coefficients are functions of the number of terms included 
in the series for the logarithm. However the coefficients converge. 
Further the first two coefficients are xmaffected. We may now construct 
a [N, N+1] Fade approximant such that Piston theory is satisfied in 
the limit s The coefficients of [N, Ntl ] Fade approximants may be 

obtained in terms of by solving a linear system. Thus one may con- 
struct a sequence of [N, IH-1] Fade approximants. These would have the 
same asymptotic behavior as equation (3 >5) for small valToes* of frequency 
and the correct limit for high frequency. 

It is now essential to establish the convergence of this sequence of 
Fade appi-oximants . 
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In Appendix D, a similar problem was considered. First the Nemann 
series solution of the problem was established. It was then shown that 
convergence of the Fade approximants is assiured if the kernel function is 
a completely continuous kernel. 

Symbollically we may denote the lifting surface problem as, 

w = Lp = [Lq + ikL^ + (ik)^!^ ... + (ik)^log ikCM^l ik ...)]p 


where and are integral operators. The solution to the steady state 
problem is 


^steady = w . 

The iterative method described can be symbolically written as, 

P= ^0^^ ’’ ^0^ • • •+(ik)\og(ik) (Mq+ ik •••)]? 

The above is a Fredholm integral equation of the second kind of the type 
considered in Appendix D- Complete continuity of the kernel function in 
the above equation can be assumed in view of the uniform numerical con- 
vergence of the kernel function method. As explained in Appendix D the 
condition of boundedness of the kernel fimction is not related to the 
condition of boundedness of the operator or complete continuity. 

Complete continuity ensures that the integral equation may be reduced 
to an algebraic system and that the solution converges uniformly. If 
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it is so, L may be tho\ight of as a matrix and 


P = 


l“^w = 


[adj(l)3 

det(L) 


w 


If the elements of L are rational functions or polynomials of a parameter 
s, then elements of are ratios of polynomials- It is from this 

fact that we are able to conclude that there exists a sequence of Fade 
approximants if the kernel function is completely continuous ■ Hence 
the theory developed in Appendix D may be applied to the above system 
and it could be concluded that there exists a seqiaence of [N,lH-l] Fade 
approximants that converge to the exact pressiire distribution as N ®- 
The same conclusions hold for the aerodynamic loads- 
Thus we may write q. .( s,M) approximately as. 


XJ 


N+1 , N 
as T a, s 
o 1 

M N , N-1 

Ms T bgS 
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In compressible flow a^ may be calculated exactly from piston theory - 
In the case of M = 0, we can easily verify from the equation that the 
loads are given by a sequence of [N ,N+2 ] Fade approximants - In the 
1 1 Tni i. M — ^ 0 and s we have the integral equation for virtual inertia 
given by 


2 z(x ,y) 
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This is a self adjoint integral equation and can be solved by a kernel 
function technique. Unlike compressible flow the pressure distribution 
has no singularities along the boundaries. It is important to realize 
that virtual inertia and Piston theory loads are independent of the 
direction of fli^t. However as circulation begins to develop the Kutta 
condition is of importance and the classical theory of motion of a body- 
through an ideal fluid is not valid. W. P. Jones [50 ] has proposed a 
simple method for solving the above equation for rectangular wings . Two 
dimensional broadside configurations such as circles and ellipse with 
symmetrically placed fins have been treated by Bryson [5I, 52, 55] with a 
view to analyzing the stabilizing effect of control surfaces on aircraft . 

For a circular wing of radius r^, in incompressible flow it could 
easily be sho-wn (Appendix C) using the same technique as van Spiegel's 
[5^] solution for low frequencies ^ that for large frequencies 


ACp(r,0) = O(k^) 

The actual computation of the Pade approx imants , though possible 
in principle, was not carried out. Exact Pade approximants could also 
be computed in compressible flow- However the amount of computations 
invol-ved to get meaningful results is enormous, as it is essential to 
expand several special functions and their products in power series . 
Also circular and elliptic wings are not used often in practice . The 
numerical technique outlined in the previous sections with some modifi- 
cations for three dimensional lifting surfaces is more useful. 

For three dimensional lifting surfaces the aerodynamic loads may 
be written as 
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q.j(s,M) = (M) + sq^^ (m) 


+ 


CqS + C^S + . . . Cjj_^ 


s + ^22 


W-1 


.. b 


W-1 


2 

s 


where q (m) -f sq.^. (m) are the exact loads in the limit s 0 and 
^ - q?j(M)-M. qj_^j (m) and q^^ (M) may be calctdated from 

and PjQ . 

If [Aq -1- ikA^] is the influence coefficient matrix for small values 

of k and {w.} = {w -f ikw }, the downwash vector, the pressure 
J J ^ J 

coefficients for small k are given by 


p . = [A ] "■" [ {w } - ik A^ [Aq ] ”^w ] 


[Aq ] may be calculated exactly as in steady flow- [A^^ ] is calculated 
from the coefficient of ik in the expansion for the kernel function 
given by equation Bl8, using either the doublet -lattice or kernel function 
technique - Integrating the pressure with appropriate displacement mode 
shapes leads to q ( m) and q. . (m) • 

The numerical evaluation of the coefficients b^ and c^ in compres- 
sible flow is similar to the two dimensional case- In incompressible 
flow since the numerical technique for the calculation of the virtual 
inertia loads for arbitrarily shaped wings was not completed, a^ was 
assumed to be unknown and bg euqal to unity - 
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After some numerical experiment N was chosen to be equal to 2. 

This leads to a first order rational function in incompressible flow. 

A simple lag is sufficient to model the circulatory effects in this 
case. In compressible flow two lags are essential^ for large aspect 
ratio wings at low mach numbers , one to account for the lag due to 
compressibility effects at high frequencies and the other to model the 
lag due to circulatory effects at low frequencies. For moderately low 
aspect ratios at high subsonic mach mmibers the lag due to compressibility 
effects may be modelled by a single lag. For low aspect ratio wings 
(AR < l) there exists no physical reason for aerodynamic lag and hence 
no lags are required at all Mach numbers. Some examples are presented. 

(1) Circular Planforms : A circular wing is considered in 

incompressible flow. For circular wings the exact results of Van Spiegel 
and Benthem and Wouters [55] for simple harmonic motion in plunging^ 
pitching^ chordwise bending and spanwise bending modes are used. 

Rational function approximations to these loads are given in Table ( 9 ) • 

The error in the virtual inertia lift in plunging is about 50^ . This 
high value is due to the lack of any data points for k > 1 . The growth 
of lift and moment for a timewise step change in the downwash in the 
plunging and the spanwise bending modes are shown in Fig\ires 22 to 25- 
No other results were available for comparisons - The corresponding indicial 
functions and the indicial functions for the lift and moment due to a time 
wise step change in the displacement in different modes are tabulated in 
Tables 10 and 12. The lift and moment in the spanwise bending mode 
reached the steady state value much faster than the chordwise bending 
mode. This is to be expected as in the first case^ the major contribution 
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to the load is from the pressure at the wing tips where it does not change 
aprreciahly with time • 

( 2 ) Elliptic Planform; An elliptic wing, AR = 6 , is considered 
in incompressible flow. R. T. Jones considered elliptic wings of aspect 
ratio 6 and 5 and computed the indicial response for the lift due to 
plunging. In the limiting case for high frequencies, he obtains approx- 
imately , 


q.]_3_(s) - tt( 




) I (s" + s(l + 0( i ))+0(J)) 


where E(*) is the complete elliptic function of the second kind and 
represents the ratio of the semi -perimenter to the span. The first 
term represents virtual inertia. It may be noted that it is independent 
of the direction of flight and hence is minimum for circular planforms . 

For low reduced frequencies, we may compute the asymptotic value 
of q^j^ by the iterative technique described, using the kernel function 
method to solve the integral equations. Fifteen spanwise integration 
points, three chordwise collocation points and ten chordwise integra- 
tion points were used for the calculation. Thus for small s we have. 


(s) = 


sq 
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+ s q 
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qj, is identically equal to the stability derivative C may be 
approximated by the formula. 


27T 


(PAR) + 1 .T:^ 72 (PAR) 

'o ^ ( 0 AR)^ + 5 . 7680 (PAR) + 6.911-89 
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The corresponding formtiLa for rectangular wings is. 
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Both the above formulae have the correct limiting solutions for AR oo 
and 0 . 

2 ^(y 

^11^” V- ) is also computed for elliptic and rectangular wings 
for increasing values of AE and shown in Figure 26. The curves indi- 
cate a maxima for circular and square lifting surfaces . For lower 

aspect ratios, low aspect ratio wing theory is a good approximation. 

1 2 

For a circular wing q and q may be obtained from the exact cal- 
culations of Van Spiegel . The errors in the numerically calciiLated 
values, using the method described above, were found to be 0.l6^ and 
0.9^ respectively. 

From these asymptotic values it is possible to construct exact 
[Ij 5] Fade approximants . Taking the Lapalce inverse transform we 
obtain the indicial response to a impulsive displacement. These results 
are compared to those R. T. Jones in Figure (27). Though the difference 
in the two results is not too large for high frequencies ( 1 ) it is con- 
siderable in steady flow. 

Next the lift due to an elliptic wing, perfonning simple harmonic 
oscillations in the plunging mode was calculated at ten frequency points 
using Laschka's collocation technique with fifteen spanwise integration 
points, eleven chordwise integration points and three chordwise colloca- 


tion points. The least squares technique for computing Fade approximants 
described earlier was used and the indicial function was computed for an 
AR = 6 lifting s\irface . This is compared to the results of R. T. Jones 
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in Figure 28. This result compares very well with the corresponding 
indicial function, obtained by using only the asymptotic values. This 
is quite significant for as it clearly indicates that with only the 
calculation of the asymptotes and the aerodynamic loads for one 
frequency of oscillation, it is possible to predict aerodynamic loads 
for all types of arbitrary motion. 

(3) Rectangular Planform: Rectangular planforms of aspect ratio 

4 and 6 , in incompressible as well as compressible flow, are considered. 
The aerodynamic loads for simple harmonic motion are obtained using 
Laschka’s method with fifteen spanwise collocation points, three chord- 
wise collocation points and eleven chordwise integration points- The 
optimum number of chordwise integration points is the integer value of 
[n(m + 0 - 5 )], n = 1 , 5, 5^ 1 , •••, m ~ ntunber of chordwise collocation 
points* It has been shown by Rowe [24] that not using the optimum 
number of integration points could lead to considerable amount of 
oscillations in the predicted values of the loads as compared to the exact 
values. The numerical technique is used to calculate Fade approximants 
and indicial functions- The indicial functions for plunging are compared 
with those of W . P. Jones [11] for AR - 6.0 and 4.0 in incompressible 
flow (figures 29 ^ 50;, and 51 Table 11). In [ 56 ] the result for AR = 6.0 
is compared to the one obtained by a potential flow technique . The fact 
that the virtual inertia for an AR = 6 wing is slightly lower than the 
case AR = 4.0 cannot be attributed to any physical reason. In fact it 
can be argued that the virtual inertia must be a minimum for a square 
lifting surface (AR = l) and have the same value as a two dimensional 
airfoil AR ~ 0. 
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When compressibility effects are included (Figure 52) the influ- 
ence of finite aspect ratio is to reduce the steady state lift while 
the starting value given by Piston theory is unaffected. Sweeping the 
surface backwards, with the planform area a constant, further 
reduces the steady state lift without effecting the starting value « 

In Figure (55) and Table 15 is shown a comparison of the 
generalized aerodynamic loads due to plunging (l) and pitching (2) for an 
aspect ratio 6 wing. This indicates that the basic behavior of the loads 
due to the pressure in any one of the modes is about the same. Thus 
this raises the possibility of approximating all the loads due to the 
pressure in any one mode, by the same combination of exponentials 
with different magnitudes . Also one may conclude that the exponential 
rise times in all the loads due to the pressure in the same mode are 
interelated. This fact led to a simplified scheme of approximating 
the aerodynamic load matrices for purposes of flutter analysis . This 
scheme is briefly discussed in the next chapter. 

In Figure (5^) the influence of compressibility on the indicial 
response in the plunging mode for a finite aspect ratio wing (ar = 6.0) 
is shown. This again is similar to its two dimensional coimterpart 
shown earlier. The steady state lift is considerably reduced due to 
induction effects . 

The indicial rolling moment for an impulsively rolling wing is 
shown (Figure 55) for a straight and swept wing of AR = 6. Here the 

influence of sweep is to reduce the lift unifoimly for all value of 
time . 

(4) Tapered Planform: A tapered planform of AR = 5.84 and a 

taper ratio of 0.524 was considered. The indicial lift for plunging 
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was compared with that of W. P. Jones in Figure 36, and indicates 
good agreement. The influence of taper is to raise the magnitude of 
the indicia! lift quite uniformly for all time. This is to be expected 
as the pressure falls to zero, near the wing tips euid hence does not 
contribute much to the total lift. Hence the influence of tapering 
with AR and planform area fixed would be to increase the total lift, 
both in steady and unsteady flow. 

(5) Swept Wing (A = 45°) with a control surface (AH = 4.l6): 

This wing (Figure 37) 'was studied extensively using the improved 
version of the doublet -lattice method described in the appendix to 
calculate the aerodynamic loads on the wing for simple harmonic 
oscillations • 

The wing was modelled as a 'beam-rod and its mode shapes (first 
bending and first torsion) were calculated by curve fitting discrete 
mode shapes (Table 14) . These two polynomial mode shapes and a rigid 
control surface mode were used to calculate the aerodynamic loads 
on the wing. 

First the convergence of the doublet lattice method in the limits 
of small and large reduced frequencies was investigated. For low 
reduced frequencies, the method was applied with l) 6 chordwise boxes 
(NCB) and 5 spanwise boxes (NSB) 2) NCB = 6, NSB = 10 and 3) NCB = 9, 
NSB = 5 (Table 15) • 

This indicated that increasing the number of spanwise boxes 
(from 5 to 7, 8 , and 10) resulted in a converging resvdt. However 
changing the nmber of chordwise boxes did not alter the results signi- 
ficantly. 
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To verify the iisefulness of the box integration scheme at high 
frequencies it was assumed that the pressure at the sending point 
(quarter chord point) of each box is given by 4 /m times the downwash 
at the receiving point (three quarter chord point) of each box. In 
this manner the generalized forces are computed in each of the three 
cases given above and compared with the exact piston theory results • 
This did indicate convergence with increasing number of chordwise 
boxes. These results are shown in Table l6. 

Finally the method was used to calculate the loads at a number of 
frequencies starting from small values of frequencies to high reduced 
frequencies (up to 5 based on semi-chord). 

The generalized forces due to the pressiire in the bending and 
control surface modes did converge slowly to the piston theoiy result. 
However the generalized forces, associated with the torsion mode did 
not converge to the resxfLts of piston theory. Hence using the results 
of the first case and the numerical method for calcxilating Fade 
approximants , it was found that the poles of the approximants were 
unstable in the loads associated with the torsion mode. In the sec- 
ond case (Table 17) the results did improve considerably at low fre- 
quencies but not at high frequencies. It was not possible to con- 
struct a stable first order Fade approximant for the generalized 
force in the second mode due to the pressure in the first. Also the 
rise times for the loads (Figures 38 - k-6) associated' with the torsion 
mode were small when compared to the first bending and control surface 
modes. The partial moment on the partial span flap has the same 
behavior as the corresponding indicial function for two dimensional 
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aerofoils . Also since the results for simple harmonic motion did 
converge to the piston theory limit in this case the indicial function 
seems qiiite accxirate • However there is no check on the accuracy of 
the results as no comparisons are available . 

From equations ( 5 - 5 ) and ( 5 -^) we may conclude that the influence 
of finite aspect ratio would be to push the poles of the approximants 
further into the left half of the 's' plane. For this reason eiLso, 
it is important to calculate the aerodynamic loads for higher reduced 
frequencies in order to get acciirate representations of the loads 
than would be necessary in the two dimensional case- An efficieint 
and fast method for calculating unsteady aerodynamic loads at al.1 
frequencies is therefore essential. 

On the other hand for aeroelastic purposes ^ it may not be necessary 
to enforce the piston theory limit for low aspect ratio wings . How- 
ever extensive flutter calculations are essential for any concliisive 
results . 

The doublet lattice method does not converge uniformly for all 
frequencies • The reason for this appears to be the lattice integra- 
tion scheme • Hence the doublet lattice method is not too useful for 
finite state modelling of aeroelastic systems . 
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SUGGESTIONS FOR FUTURE RESEARCH 


Generally for purposes of flutter analysis, it is sufficient 
to calculate aerodynamic loads accurately at the frequency and 
flight velocity at which flutter begins to occur. However for 
designing control systems for the suppression of flutter, it is 
not only essential to predict the root locus in the neighborhood 
of the flutter frequency, with respect to the flight velocity, 
but also some of the stable roots of the transcendental characteristic 
equation which influence the unstable locus. Thus aerodynamic 
models used for flutter suppression studies must be fairly accurate 
over a range of frequencies. Numerical techniques must converge 
uniformly over a large bandwidth. Collocation techniques discussed 
earlier for the calculation of aerodynamic forces generally satisfy 
these convergence requirements. 

The aerodynamic modeling theory presented in this report 
suggests the possibility of representing the aerodynamic load matrices, 
Q^, in a certain modal coordinate system X, as 


Cxij* ■ 


{n . . (s) /r . 


(s)} = 


P R 

X X 


-1 


where R is diagonal. Transforming to the actual modal coordinate 
system Z = TX 

Q = t"^Q^T = t"^P^T [t"^R^T ]"^ = P(s)R(s)"} 
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where R(s) is not diagonal. Hence in general 


P(s) = E V"" %+l 

n^O 

and 

N-1 

r(s) = E R s" + MIs^ 
n=0 

The interpretation of this approximation is based on the fact that the 

loads due to the pressure in any mode may be approximated by the same 

combination of exponentials with different magnitudes. If M = 0, 

this reduces to a [N + 1, N - 1] matrix Fade approximant. The 

coefficients R^ may be obtained numerically using a combination of 

Lagrangian interpolation and least sqtiares technique • This has already 

been verified for two-dimensional aerofoils and in certain special cases in 

three dimensional wings . In these cases , however, the Piston theory limit was 

not enforced and N = 0 and N = 1 were sufficient in the range of frequencies for 

which the phenomenon of flutter becomes a possibility. 

The above approximation of the aerodynamic load matrix has 

several advantages over the method of approximation already discussed. 

The matrix approximants lead to only (M^(N-I) in compressible 

flow) states, in addition to states required to define the structural 

system, where is the n\miber of vibration modes considered for 

modelling the structure . Approximating every individual element of 

2 2 

the aerodynamic matrix by a Padd approximant leads to ‘N (M^ (N-1) in 
incompressible flow) additional states . 

The equations of motion of the entire system can be written as 
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- 1 , 


- P^P(s,M)r" (s,M)qj^ = BQ^ 


( 1 ^. 2 ) 


\ - H^q^ + H^q^ 


P„ = 


q A b 
^ w 


The equation reduces to, 


[M. 


[j^RCs^M) + K^R(s,M)- P^P(s,M)]R"\s,M)q 3 ^ = 




The flutter s^eed is determined hy 


2^2 

det[M^R(s,M) + K^R(b,M)- P^P(s,M)] = 0 

b 


(4.4) 


The eqiiations of motion can be reduced to the form 

(i^-5a) 
(4.5h) 


X = FX + GU 
RX 


where X is the vector of states {y]{y} etc. and U are the control 
torques on the aerodynamic control surfaces . There are several 
criteria when such a realization would lead to a system iidiich is con- 
trollable and observable [57 ] • These are generalizations of the 
criteria for scalar transfer functions . The survey article by 
Silverman [58] presents an excellent description of the realization 
problem . 
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The model is useful for synthesizing control laws for active 
flutter suppression. Control laws may be written in terms of the 
state vector [x] as 


U = - CX 

/N. 

C may be obtained by arbitrary pole placement. The characteris- 
tic polyncaniaJ. of the open loop system is assumed to be 


X(f) = s“ + Z a.s“~^ 
i=l 


and that of the desired closed loop system. 


X(F - GC) = s" + Z 

i=l ^ 


The coefficients a^ are related to F by the relations, 

1 ^ 

\ = k .fo “j 'v+l-j ’ i = 1. 2, 3, ... N 
'^ere = trace (F ) and a^ = 1- Similar relationships hold for 

✓s, 

(F - GC) and p. . Using the notation, 


T A r T 

P ~ ^ y • • • P^ J 


T A 

a = [a^ ^2 ^ ^ 


and 
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A i 


1 


0 

1 


0 

0 


0 

0 


0 

0 


a , a ^ 
n-1 n-2 


n 


* n 

n-1 


1 


0 

1 


it is possible to obtain a relationship between the control gains C and 
[p]-[a]. It may be noted that A"^ always exists and can be evaluated 
easily* In [59] it is shown that provided F is cyclic (if not an 
initial gain may be chosen such that [F - GC^] is cyclic] 


trace GC 

/N 

trace FGC 


_-l 

= A (p-a) 


trace F 


n-1 


GC 


This is a linear system of equations in terms of the unknowns in C and 
may be easily solved* Preliminary resiolts are extremely encotxraging- 
A detailed discussion of the modeling procedure and numerical results 
is beyond the sco]^e of this paper, A simple scheme for control 
law synthesis is shown in Figure ^7* This technique is outlined by 

Lyons ^ et al* [6o]. Optimal laws based on quadratic synthesis may be 
obtained by the eigenvector decomposition technique as applied in 
Ref. 61. An offshoot of this technique would be a new method 
of flutter analysis which would be useful In structural optimization [62]. 
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CONCLUDING REMARKS 


A general theory of finite state modeling of the aerodynamic 
loads on oscillating airfoils and lifting sxirfaces has been systema- 
tically developed and presented. In particular, the natiare of the 
behavior of the unsteady aerodynamic loads in the frequency domain 
is now well understood. The analytical reasons for the difference in 
the behavior of these loads at high frequencies in incompressible and 
compressible flows are explained. Extensive studies need to be 
carried out for supersonic lifting surfaces » where it is expected that 
higher order Fade approximants would be necessary for an accurate 
prediction of the aerodynamic loads for arbitrary motion. The theory 
presented has several useful applications, especially in the control of 
the vibration modes of aircraft wing structures and in the active 
suppression of aeroelastic instabilities. 

This study sioggests the possibility of approximating aerodynamic 
loads as Fade approximants in Mach number as well. However, in the 
transonic regime nonlinearities need to be considered to account for 
shock waves . 

Several other suggestions for improving the results and for 
applying the theory have also been presented. 
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APPENDIX A 


INFLUENCE COEFFICIENTS FOR A BEAM-ROD WITH A FLEXIBLE ROOT 

The inflixence coefficients for bending of a beam of length Z are 
given by (the deflection at point x due to a unit load at point 0 

^ > X > 0 
£ > X >§ 

where EI(a) is the flexural rigidity distribution, Is the root restraint 
stiffness constant for rotation in the plane of bending, and is 
the root stiffness constant in shear. The influence coefficients for 
torsion of a rod of length Z are given by (torsional deflection at 
point X due to unit moment at point $ ) 

' I ?>x>o 

0 0 

>diere GJ(X) is the torsional rigidity distribution and K is the 
torsional spring constant at the root. If the root is rigidly fixed, 
the influence coefficients reduce to those given in [l]. 


(^-X)(x-X) , ^1 

EI(X) ■ 


^ (§-X)(x-X) , §x 1 

= J iJrTTV" » 
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APPENDIX B 


INTEGRAL EQUATIONS FOR UNSTEADY AERODYNAMICS 


Three Dimensional Lifting Surface 

The source solution of the linearized equation for the velocity- 
potential is given by 

4-tt^ = ~ (B-1) 

where H(-) is the Heaviside step function, r is the distance between 
the observation and the source point, x is the time measured from the 
instant of the disturbance and c is the speed of sound. The observa- 
tion point is affected at time x if the disturbance moving with a 
velocity c can traverse the distance r in that time; the strength of 
the disturbance is -l/hirr. 

A coordinate system (C, 11 , Q, t) moving with respect to a fixed 

coordinate system (xQ,yQ,ZQ,tQ) , along the negative x^ axis with a 

velocity v, fixed in the wing, is considered. The fixed and moving 

coordinate systems are assumed to coincide at time t = t. The 

o 

coordinates of the point ( ^, Tj, Q, t) referred to the fixed axes are 
(€ - ^ ^ acceleration potential source at the 

point (C, 'H, £ ) in the moving reference frame will produce an 
acceleration disturbance at the point, (x, y, z), after some time delay; 
the magnitude of the disturbance is 
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2 

= (x - C + v(t^- t))^+(y - Tl)^ + (z - C)^ (B-2) 

In fixed coordinates, the velocity potential is related to the accelera- 
tion potential as 

tg _ _ 

^(x,y,z,t) = J* 9^(x,y,z,t)dt 
^1 

idiere [t^^t^] is the time domain of influence of ^ich affect^ the 
velocity potential at time t in fixed coordinates, t^ and t^ are the 
last and first instants of time at which the moving source can affect 
the observation point- 

Therefore the velocity potential at (x,y,z,t) due to a moving 
doublet at (^ , T], Q, t^) of strength t^) is given by, 

§(x,y,z,t) = - I T1, C, t^) I dt^ (B-5) 

The time satisfies the equation, 

((x - ^) + v(t2~ t))^+(y-Tl)^+(z-C)^ = c^i^ 2 ~ 

This is the equation for en ellipse and determines the region of 
integration for equation (B-5)- IfX=x-6,Y = y- T) and Z = z - £ 
and t^ - t = T, we have 
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MX + R 


1 - M 


, R - MX 

+ ^ 

M - 1 


where R + (l - M^)(Y^ + Z^) 


For M < 1, only one of the solutions is negative- Hence 


^ ^ 1 MX - R 

^2 “ ^ c , ? 

1 - M 


If the perturbations in the flow field existed for infinite time in 
the past then t^ = - <“ . On the other hand if the flow field is 
tmperturbed for time t < 0, then clearly 


^ ^ L 1 MX - R 


t^ = 0 provided 


1 - M 


1 MX “ R 

If t < ~ -R y t, = - 0^ as the disturbance has not reached 

1 - M 

the observation point. For M > 1;, both solutions for T are negative 
and the period of influence is the time an advancing spherical wave 
woxild take to pass the point of observation. Therefore^ 

t _ t + i R - MX ] 

2 ^ ^ M^ - 1 


^ ^ 1 R + MX 
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provided, of coTorse , both and are positive. For M = 1, (transonic 
limit), we may let M=l + eorl-c and take the limit as t — ♦ 0. If 
we denote = X + V(t^ - t), i = 0, 1, 2, then in subsonic flow. 




, = X - Vt 


provided the flow field is tanperturbed for time t < 0 . 
For supersonic flow, t s\rfficlently large. 


(B-9) 


„ MR - X „ (MR + X) 

Ug - -2 ^ 2 

M - 1 M - 1 


(B.IO) 


Hence the velocity potential at (x, y, z, t) due to a moving 
doublet at ($,”n,C,tQ) is equal to 


§(x,y,z,t) = - 


1 a |- 

■577V ^ J 


II X 

^2 + 


V + 


dU() (B.ll) 


The total velocity potential at (x, y,z,t) due to a distribution of 
moving doublets on the lifting surface S, at 5 ^ 0, is given by 

II X 

4-7r$(x,y,z,t) = - | J‘J‘ d? dn ^ J‘ ^ dU (B. 12 ) 


“1 ^0 


2 2 2 
+ Y + Z 


The downwash on the surface Is given by 


. a u A(^(5,n, o,t + u -x)du 

4. w(x.y.t) = - i ^ dC dy A ;2 jr 

^ </ 2 2 T 

U + Y + Z 

O 


(B.13) 
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It can be shown that a discontinuity across a doublet layer exists 
such that fl] 


A(|>(x,y,0,t) = $(x,y,0,t) - 'l>(x,y,-0,t) = 

(B.14) 

where Ap is the pressure differential across the surface. 

Substituting (B.14) in (B.13) and performing the differentiation 
using Liebnitz rule, we have the following integral equation for the 
downwash in subsonic flow. 


MX+R, 


0 ' ^ UCT-T) 


„ dU 
. 3 o 


U +Y 
o 


+ / Ap(5,n, (t-T)) 

0 ((X 


H(t-X/V) Vdxl 

X-Vt) +Y ) 


X-MR, 


where Ap(5.n,t) = 0 for t<0, U = , 

1-M 


2 n . ^ 1 

q = Jsp V^,R, = hr + (l-M^)Y^ and t > ~ • 
“ 1 V c 1_M^ 


(B.15) 


2 2 

In the limit as t 0, U -► X, MX ^ R^ and X + Y -> 0, and only the 
first of the 3 terms remain in the equation B.15. Integrating over a 
small region around the receiving point (x,y) , we have 


w(x,y ,0) _ _ M Ap(x,y,0) 

V 4 q 
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which is the familiar piston theory resiilt. At the initial instant, the 
lifting surface behaves like a piston, and the flow over it is 
similar to one dimensional compressible flow. It may be noted that a 
similar equation obtained by Drischler [65] may be reduced to the 
same form as (B-15) after some algebraic manipijlations . 

Taking the Laplace transforms of Equation B.I5, we may obtain an 
integral equation relating the Laplace transform of the downwash 
w(x,y, and the pressiire Ap(x,y, ^) , where p is the Laplace trans- 
form variable • 


btt 

V 


where K^(x,Y, ^ ) = 


q de dT] Ap(c,Tl, ^)k^(x,Y, ^ ) 


^ X - MRt 


M 


MX + 




\ . „2 


X + Y 


X-MR, 




0 


(U, 


2 2372 + ^ 5/2 ^ 1 

y^)^/ 2 0 ((x _ vt)^+ y2)5/2 J 


(m) represents the Laplace treinsform of m. The third term converges 

only for real (s) >0, and ^ 

P5"r . p ( ) MX + R 0 

^ V ^ 1 - m 2'' ^ ^1 . p ^ 

% X^ + Y^ MR^-X (t^+Y^)^ 

1-M^ 


K^(X,Y,^) = e "" 


PI 

V 


dT 


JL ( E!i1) 

„2 '' 2V '' 




P Y 
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■where |hI j_(‘) is a Struve function and 

Y(‘) is a Bessel function of the second kind. 

For simple harmonic motion this reduces to the kernel function given 
hy Watkins, et al • . For purposes of calculations however it is 

convenient to write, Kj(.X, Y, ^ ) as. 


Kj(x,Y,^) 


V 


^ MR -X 

1 -m'*^ 


M ^ + % 
^1 X^ + Y^ 


00 


+ J 

MR^-X 


1 


-pT 


e^ dT 

(t2+ y")5/2 


As shown hy Landahl [ 65 ], the integral is also equal to. 


(B.I 7 ) 


Y MY -XE^ 

~2 2 
X + Y 


-PYt 




dr . 


However equation B.17 is most convenient for calculation purposes- The 
second term in equation B-17 is a singular integral at Y = 0- Hence 
the integration over the lifting surface must be done in the sense 
of Hadamard, as expalined by Mangier [66]. This is indicated in 
equation B-l6 by the symbol 
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For small values of the reduced Laplace transform variable s = pb|v, 
we may write the Kernel f\mction as [6j]. 


K^(X,Y ,s) - 2 (a s + s ^n(s)"b s^) 

^ it=0 "" ' ' n ^ 


(B.18) 


where i (i + E ) 


1 1 


2 ? 
X + Y 


\ , 


^0 = V2 


^ 2b 


The coefficients b^ are functions of the spanwise variable Y only (and 
not of X) 

In steady flow, for s = 0, 


K,(X,Y,0) = (1 + 

If VA 

and in the limit as s from piston theory 


) 


(B.19) 


= 2t7M6(X,Y) 


(B.20) 


In the limit as M 0, we have. 
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K,(X,Y,s)L_^ = J' 


- (T+X) 

00 ^ b 

e 


^ -X (T^+ y^)V2 


dT 


(B.21) 


In two dimensional incompressible flow we have ^ 


K^{X,s) 


= K_(X,s) = J‘ J 


-s(t+X) 


M=0 20 


X (t^+ 


dT dY 


_ 2 2s -s x/b 


X + ^ e “ E.(sX/b) 


= I - ^U(1A.- ^ ) (B.22) 

where is the exponential integral and U(l^l^-sX/b) is a Kumraer 

function. 

The above Kernel function indicates a distribution of time delays 
over an infinite time domain. Further they are not valid for Real(s) 

< 0. However if a solution to the integral equation can be obtained 
for Real (s) >0^ the convolution principle may be used to find the 
solution for any s . 

For purposes of calciilations , we note that the function 

|h| ^ (.5-1^ ^ ) in equation B-17 may be approximated by a 

convergent sequence of rational functions in or Fade approxi- 

mants . It follows therefore that — — 5-7— can be approximated by a 

(r^+1)'^/^ 

finite sum of exponential functions . Two such approximations are 
known, one due to Laschka [25] and the other due to Runyan and Watkins 
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[68]. The errors in these approximations are compared in [6 q]. 

ar 


Laplace transforms of exponentials of the type e are given by 
rational functions* Hence ^ in principle^ it is possible to approximate 
the entire Kernel function by a rational function in s- Further in 
order that it converge to the piston theory kernel for large S, the 
approximation must be [N, N] Fade approximant. This point of view 
has useful consequences , regarding the behavior of the aerodynamic 
loads in the frequency domain. 

In supersonic flow the kernel function may be easily shown to be 
r sU ^ X+MR^ 


Kj(x,Y,s) = e 


sX 

b 






■1 e 


\ X^+Y^ 


M 


X-] 


■MR^ 


b m2_i 


M^-2 


X-MR 


= 0 


for X > 1 |y I 


for X < (y I 


(B.25) 


In this case the kernel function is valid for all s, including the 
left half of the 's' plane. 

In transonic flow though the linearized theory' cannot predict 
shock waves . We may obtain the kernel fvinction by proceeding to the 
limit as M 1 from the supersonic case. We have, 
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= e 


sX 

b 


-X 


£U 

b 


dU 




2X 


_ ^ oi?+?)^ x^+ 
2 ^ 21 . 


X > 0 


= 0 for X < 0 (B.24) 

In this case again the kernel f\mction is only valid for real (s) >0. 

Two Dimensional Airfoils 

In two dimensional flow, the integral equation relating the down- 
wash and the pressure distribution may be obtained by integrating the 
corresponding three dimensional kernel function with respect to the 
spanwise coordinate from - ® to + For the case of arbitrary motion 
of an airfoil , one can show after some computation, that 


877w(x,t ) _ 
V 




2de 


x+T 

J Ap(C.t- 

-X 

L 1+M 


dr 




2 

M T 

(•H-X)^ 


-X 

1+M 




-X 

1-M 




(rfX)^ 


X+T 

V 


)) 


(B.25) 
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vhere = x - 


Vt(l + m) 


M 






if X - Vt iaii < Xj_ 


X _ V + - M) 

2 " M ^ 


if X + Vtii^ < x^ 


= X„ 


if X + 


Eq\iatlon B.25 is valid for any type of airfoil motion. In incompressible 
flow we have , 


8irw(x,t) 

V 


, Vt-X 

- I ^ 2 J* Ap(C, t - 

S, 


(x+t) 

V 



At time t = 0, we have 

w(x,0) M ^p(x,0) 
V " ■ ^ q 


(B.26) 


(B.27) 


Non-dimens ionalizing the time variable t to t = Vt/b and taking the 
Laplace transforms of the integral equations , we have 

"v"" ~ i Kg(x - e, S) ■^(5,s)d^ (B.28) 

where w(x,s) is the Laplace transform of the downwash and £p(x,s). 
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the Laplace transform of the pressure 


^ I V „ *nf +s*mSc 
KgCx^s) = 2s Pie 


[(MK^(Ms* 1 x1)+Kq(Ms* |x|) 


1 * 


+ e"® ^P lTi{ i^)+ sX ]■ e® K^(Ms* ]x jX)dX]3 (B-29). 


*)(■ 2 

where s = s/P , Kq(-) and K^(‘) being modified Bessel functions of the 
second kind. 

Also K^iX.s)! + 2 Mtt6(X) (B.JO) 


and 


KgCX^s) 



0 “ 


2p 

X 


(B.51) 


In incompressible flow, equation B-29 reduces to equation B.22. Equation 
B.22 has a well known inverse given by, (after non-dimenslonalizing 
all distances by the semi -chord). 


Ap^(x,s) = r G(x,z,s)w(z,s)dz 


(B.52) 


where G(x,z,s) = " [sA(xj 




A(x,z) = ^ log 


1 T 1 - xz + A - X - z^ 

1 - xz - A -y^ A 


(B04) 


Kj^(s) 


' Kj^(s)+ Kp(s) ' ^ ° 


(B.55) 

continued 
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2sU(l|, 5, 2s) 

2sU(l|, 3, 2s)+ U(i, 1, 2s) ' 


real (s) > 0 


(B.35) 

continued 


= s [1 


- J {(Kq(x)- K^(x))^+ /(Iq(x)+ I^(x))^} 


-xt 

e 


2 

X 


dx] 


K^(-) and Ij(') are modified Bessel functions and U(a,b,z) is a 
Kummer function. 

Equation (B.28) may be reduced to several alternate forms. One 
method of doing this reduction was given Fettis [70]. The integral 
equation may be reduced to a non-singular Fredholm integral equation 
with a ccanpletely continuous kernel function (defined in Appendix d) 
given by. 


^ (x,s) * p 1 

Ap(x,s) = p s M J'^K^^(x,C,s)^(€,s)d€ 

L.(x,e,s) = H(c - x)- ^2|ii . 

TT TT V 1+x 


(b- 36 ) 


16 ”*' 

” ~2 ^ ^(^^z,s)k(s (z-C))dz. 

SM -1 


H(€ - x) is the unit step function 


and 


K(s(z-{) = K,(z - ?,s)- . 2 


■ 


‘2'“ p 

Another important representation of the integral equation may be 
obtained. The kernel function Kg(X,s) may be written as. 
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K^(X,s) = 3{K2 q(X,s )- ^^[K^^(X,s )-K^^(X,s )] 


2 " 20 


20 ' 


+ !^.3*.ti|l„ s*- 4 “4 in I - 2V1 

M 


-•s 




+ 0(fr)-0(X^+ aX log X)} 


(B.57) 


= e'*'^ ^ P(K2 q(X,s*)- ^[K2q(X,s*)-K2®q(X,s*)] 


where 


+ (l-2[ln(|Ms*)+ y])+ M^K^(X,s*)} 


Kg^X^s*) = 2(| + s*ln X) 


(B.58) 

(B.59) 


K ^{X,s*) = 0(M° + a In M) 


(B.UO) 


Using equation B-58 we may invert eqiiation B-28 using the well known 
solution for incompressible flow. Thus we have 


q Y 27T 

1=0 




+ B^(st x)j^(s*,x)]- ^ e^® * J J' G(x,z,s")K^( 


J J G(x,z,s*) 1C (z-C,s)e“^® 


- 1-1 


dCdz 


(B.hl) 


^diere 
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dz 


(B.te) 


AP(^(x,s*) 


= f J G(x>Z;S*)w(z,s)e"^ 
^ -1 



jV ac . 

X ^ 


Aq(s ,x) = In" x(2s -s x)+ ^ (s''0(s'')ln 2-(^(s'')-l) 


■Jf 

+ ^(0(s*)+s*(l+x))+s*^(l+ In 2)(1 + x)). 


* , * 


Bq(s ,x) = -(2 - s x)tts 


A^(s*,x) = s*^ln"^x + 


/ * N *2 

B^(s ,x) = -7TS 


a = [1 . 2(ln(| Ms*)+ y) ] 


Thus for sTifficiently small M we have for Ap, 
Ap(x,s) ^ M^s*x 

q q ® 


2 2 * 

M Ms X I 

27T ^ 


2 A^(s jx)j (s,-l) + B. (s ,x)j.^(s ,x)] 
i=0 ^ ^ 


* _ s.O, * 


4 , 

+ 0(m s (a + b In M) ) 


(B.45) 
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This solution is exact to 0(M^s*). The first term is identical 
to Osborne's thin airfoil theory [71]- The usefulness of this 
representation will be briefly discussed in Appendix C- 



APPENDIX C 


CALCULATION OF AERODYNAMIC LOADS 

FOR SIMPLE HARMONICALLY OSCILLATING AIRFOILS AND LIFTING SURFACES 


In this appendix some of the existing methods and their extensions 
for the calculation of aerodynamic loads on wings with and without 
control surfaces will be briefly discussed. 


Two Dimensional Airfoils 

Analytical methods > - In incompressible flow Theodorsen’s 

solution is well known for calciiLating pressure distributions and air- 

w n 

loads on two dimensional aerofoils • For a downwash = W^x , 
we have , 


q 




1-x' 

1+x 


Wj^[x’^+( 


n-1 




£=0 


■W 


X ik(l+x) . 
'' n+1 


( 


n 

E 




- l+W'')-!) -iM 1 


, X /H ,x for n even 

where I (n) = ! ) 


nil , n-1 , 

2 ■ 2 ‘ 


for n odd. 


(c.i) 


W / 

For a downwash given by — = 

= 0 



for X > X 

c 


for X < X 

c 
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1-x 

1+x 


w 

n 


n-1 n-r-1 , 





cos^x + 
c 


ik- 


(l+x) 

tri-1 


r=0 * 


n-r 
( 2 (- 
IfcO 




n 


+ (l - C(k ) ) S C (-X ) (J + J ) 

^ r' c n-r n-r+l"^ 

r=U 


, ik(x-x ) 

W i(x-x) + V - ■■ ■■ 

TT n c' n + 1 


n+1 


}log 


1-x X + A-x ^ A-x^ 
c c 


(C.2) 




1-1 


where J 


2,g 


2?, ^ c ^ ^ 


(2|-l)(2^-5). ..(2«-2R+l) 2«-2k-l ^ 

( 7 , ) 


k=l 2 ' (i-1) (?,-2) . . . (?,-R) 


+ 


2 £ ! 

PP 

2 h H : 


-1 

cos X 

c 


and 


2i+l 




2i + i 


'x^^+ 
I c 


H 


-i k.+l 

E T • • • (i -k) 

(2^ -1 ) { 2 ^ -5) . . . (2^ -2kH) c 


2i-2k-2 


k=0 


The expressions for the loads may be easily obtained by integrating 
the corresponding moments of the pressure distributions. 

Equation B.Uj extends Theodorsen's solution to subsonic flow for 
low mach numbers and frequency. The first term of eqixation B. ii5 is 
identical to Osborne^s theory. Kemp [72] has provided the closed form 
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expressions for the lift and moment based on Osbourne's theory. How- 
ever it was found that this approach is not too attractive for ptirposes 
of n\jmerically computing pressure distributions or aerodynamic loads- 
Further these solutions are not valid for either high subsonic mach 
numbers or high reduced frequencies, when compressibility effects are 
not smsdl relative to circulatory effects. 

For high frequencies, and Mach nimibers M / 0, the method of 
acoustic planforms of Lomax, et al. [12] is viseful • The two dimen- 
sional unsteady problem for compressible flow can be fonnulated in a 
coordinate system x',y',z',t' fixed in the undisturbed fluid, by the 
equation. 


CR 


X X 


+ «P, 


z 'z 


2 'Pft' 


with the appropriate boundaiy conditions . For the three dimensional 
steady supersonic case , in body fixed coordinates , it is 


cp + cp = (M - l)cp 

yy zz ^^xx 


Hence it is seen that the solution of the unsteady problem can be 
obtained from the solution for the steady flow about a swept forward 
wing tip placed in a supersonic flow with Mach number M - v^. The 
angle of sweep is determined by the Mach nvanber of the two dimensional 
aerofoil from than relation A = taiT^/^. The wake vortices in the 
\xasteady two dimensional problem correspond to the tip vortices in the 
steady three dimensional problem. The three dimensional steady flow 
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problem may be solved by Eward’s method [1 ] • For example the pressure 
distribution obtained this way, corresponding to an indicial down- 
wash Wq, for t > 0, is given by, 

4Cp(T,x) - - g»o ' “=P1 


for T < min( 


M M 
1+M ^ 1-M 


(1-x)) 


A Cp(T,x) 


A Cp^ + 


8 \ 
M 7T 



1 /1+m)t 

\teirry 





A Cp(T,x) 


ACp + i tan"^ 

V\ M 7T 






- 1 


Vtt 

where t = -=— and 0 < M < 1 

D — 

Expanding in the variables (l-x) and (l+x) indicates that for t > 0, 
at the leading edge (x = -l) 


iC (Tj,x) « as 

and at the trailing edge (x = l) 

ACp(T,x) -»• 0 as \/l-x 
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The indicial aerodynamic load coefficients may be obtained by 
integrating the pressures over the entire chord of the aerofoil- 

Taking the Laplace transforms, and using the asymptotic expansions 
of the resulting Frenel integrals 


R -ih 

Fr(R) = J ^ du 

0 \/ 2 ^ 


1-i 

2 


+ 







) 


+ 0 (- 
R 


1 

7/2 


) 


we may obtain the expressions for the non-dimensional lift and moment 
coefficients in plunging and pitching for large values of k . 

These are 


2b -q 



^(1-M)- 


1-i 





e 


- k2M 
" 1-M 


! 


A. 

2b^q 


1-i /T 


+ 0(l/k) 


-ik2M . 2M 

e - d-M)"e" 


+ o(iA) 



■ k2M .k2M 


+ 0{l/k) 


M 

a 

2~ 
2b q 


4 ^ 

5 M 


+ , 

2 

M tjM' 






It is interesting to note that it is not possible to expand these 
coefficients in a power series in (l/k) . Also they are not valid for 
M = 0. Thus it is possible to analytically obtain the behavior of the 
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aerodynamic loads for low and high frequencies. While they contribute 
to a better understanding of the behavior of the pressure and aero- 
dynamic loads at high and low frequencies, they do not serve any other 
useful purpose. For this reason and for computational efficiency, it 
is convenient to resort to numerical techniques which are valid for 
all valxies of freq'uency and mach number. 

Numerical methods • — One of the earliest methods to calculate 

aerodynamic loads on two dimensional aerofoils is due to Dietz 
(described in Refs. 1 , 2 , and 3). For wings with trailing or leading 
edge flaps , Fettis suggests a method that is fairly accurate . From 
a computational point of view Hsu's method [73] for wings without 
flaps seems to be simplest. However Hsu's method can be extended for 
wings with control surfaces. This modification is briefly described 
below. A pressure loading function of the form, 

M o A N 

ACp= E 


log 


1-xx + 
c 




1-xx - \A-X^ 

r» “ » r* 


is assumed, where h^(x) are identical to Hsu's chordwise loading 
functions in compressible flow and are the exact loading functions 
corresponding to a downwash of x"^ in incompressible flow, given by 
equation C.l. S is the angular deflection of the trailing edge flap. 
The above loading function is substituted in the integral equation and 
integrations of those terms that contribute to discontinuities in the 
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downwash and its chordwise derivatives are performed analytically. The 
coefficients a^ can then he foxind hy comparing like terms on both 
sides of the integral equations- The kernel functions K2(x,s) may be 
written as 

--y ln|x|-(| )^(p-5P^)x ln!x| 

OO 00 

+ 2 Q + Inlxl 2 B x^. 

n=0 " rt=2 

S.Q, and a^ depend only on the terms outside the summation signs and 
are given by. 


^0 = \ ^2 = “I 

P P 

and 

a = 0( )5-0 (m) 

^ P 

In incompressible flow a^, aj^, •••, etc-, are equal to zero (from 
equation C -2) - 

The singular part of the above kernel function and the pressure 
loading functions are integrated exactly- A computer program is 
developed that evaluated the Possio kernel by Gaiissian and Berthod 
Zabrowski quadrature, integrates the assumed pressure loading functj.ons 
and the non-singular portion of the kernel function numerically and 
then solves for - The program then obtains the corresponding 
generalized aerodynamic forces for any mode shapes- Results indicate 
that for M = 0-0 and 0-7 , the error in using a ten point collocation 
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sch6Hi6 and not including the flap singularity is 20 - 25 ^ of the four 
point collocation solution with the singularity included* In incom- 
pressible flow using the four point collocation scheme results in an 
error of less than 1^, in the moment about the flap hinge line- The 
progrsm may be used for values of 0 k < 10 and can include trailing 
and leading edge flaps . For high frequencies a large number of colloca- 
tion and integration points are required. 

In supersonic flow, the expression for the pressure and airloads 
given in Ref. 1 may be integrated using Gaussian integration with no 
difficulty as control surfaces do not contribute to any singularities . 
Results for a rigid airfoil with trailing and leading edge flaps are 
given in [Ik] and [k5]- 

Three Dimensional Lifting Surfaces 

Analytical techniques . - There are no known analytical tech- 
niques available for arbitrarily shaped three dimensional lifting 
surfaces. However, as pointed out by Miles [7k], the solution of the 
partial differential equations for the velocity potential by separating 
the variables is possible only in eleven Euclidean coordinate systems- 
Of these only three define lifting surfaces of practical interest. 

These are (l) elliptic cylindrical (two dimensional eierofoil ) (2) 
oblate spheroidal (circular planform) and (5) ellipsoidal (elliptic 
plauform) . 

The circular oscillating lifting surface in incompressible flow 
was considered by van Spiegel [5k]- It is interesting to note that from 
van Spiegel's general solution, we may extract the solution for the 
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virtual inertia loading. The orthogonal coordinate system for the 
circular wing is given by 



z = uTi where 0 <t]^- 1<^<1, 0<v<27T so that the entire space 
is covered just once- The wing itself is given by t] = 0, while the 
part of the x, y - plane outside the surface is given by |i = 0. 

The problem is treated by van Spiegel by the method of the 
acceleration potential . The regular solution is equal to 

00 i -1 . . . j 

= S S (u)Q^^ (iT)){c^^cos i v + S^*^sin i v] 

i=l j^O 

where and are Legendre's associated functions of the first and 
second kinds- Since should be an odd function of z and u, the 
summation sign is restricted to i + j = odd- This is indicated by 
the prime added to the summation over j in the equation for 0. The 
constants and are determined by the prescribed normal derivatives 
at the airfoil. By expanding 0^^ = 0^^/u in a series of surface 
harmonics the constants are determined- It is foimd that^ 

1 2tt 

0^(Ti,a, ) = 0^2(d3^,Vj_)G(Tl,U,V,O,Uj_,Vj_)dV^dy 


where the Green's function of the second kind is equal to. 


G(n..,v, 0, u,,V ‘ f 

n^l m=0 m 

• Cos m( V-V ) 
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where Iq= 2, = . . . = i . 

The above solution is one that has the exact prescribed acceleration 
on the aerofoil . Further the solution is not singular at the leading 
and trailing edges and behaves as \/r where r is the distance of the 
leading edge or trailing edge. The solution that has the exact normal 
velocity w(u,v), prescribed on the airfoil Is 

* ;n ^ 27T 1 

f (l.U.v) = ( ^ + u J J^d^w(u^,v^) • 

■G(Tj,u,v;0,u^,V^)dv, du . 

This solution leads to apressure singularity at the trailing and lead- 
ing edges due to the differentation in x- 
Consider now 

0 , 0 , 

This is a solution of Laplace’s equation having zero normal velocity 
on the aerofoil and a singularity at the point r| = 0, u = 0,V =Vj^ 
which lies either at leading or at the trailing edge. The singularity 
is of the type ^(cos Q^)/r^ , where r is the distance from a point on 
the aerofoil to the singular point and a the angle between the vector v 
and the radius of the circle. Thus the singular solution is ^ 


^^cos ( V - )+ 1 ^ 


92 



3Trf2 

tp(il^u,v) = J a(v,)G (ti,u,v;Q;0,v, )dv, • 

vdiere the integration is along the leading edge only, since no singu- 
larities are admitted at the trailing edge- As in the two dimensional 
case the singularity in is of the type \/l/r, where r is the dis- 
tance from the leading edge. The function a (r^) should be deter- 
mined from the condition that the normal velocity at all points on 
the aerofoil due to the total acceleration potential agree 

with the prescribed normal velocity. Hence the normal velocity at 
all points due to the acceleration potential ^2. ^ should be 

equal to zero. This leads to an integral equation which may be 
solved by Fourier series expansions. 

On examining this solution carefully^ we find that very large 

jf. 

frequencies of oscillation or at the starting instant - ^ ) 

is independent of frequency. Hence we conclude that the solution 
corresponding to the virtual inertia loads is. 



J 

0 


s\( 




zCU;L- 


)d Vjdu 


Like the regular solution this is not singular at both the lead- 
ing and trailing edges . 

The effect of compressibility may be taken into account by means 
of spheroidal wave functions . Also the oscillating elliptic wing can 
be treated when Lame functions are used. However, these analytical 
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solutions are not really exact and the degree of computations is 
enormous when compared to numerical methods . Hence these techniques 
are not too useful. 

Numerical techniques . - Numerical techniques may be broadly 
classified into two groups (a) finite element methods and (h) kernel 
function techniques. The doublet lattice method of Albano and Hodden 
[ 26 ] and potential flow methods [I 8 ] in subsonic flow and the Mach 
box method in supersonic flow [ 75 ^ 7 ^] are examples of the first 
type. The kernel function methods of Laschka [25], Cunningham [25] 
and Rowe [24] in subsonic flow and Muller's method [ 67 ] in supersonic 
flow are examples of the latter type. Two of the methods that were 
actually used in a slightly modified form during the course of this 
work are briefly described- 

The doublet lattice method: This method for lifting sur- 

faces is based on the developments of Albano and Hodden. The surfaces 
are divided into trapezoidal panels, the downwash being taken at the 
midspan of the 5/4 chord- A distribution of acceleration potential 
doublets whose strength is to be determined^ is assumed along the 
quarter chord of each panel- The assumed pressure loading is 
then integrated along the quarter chord line- In the original method, 
the integrals were obtained by approximating the kernel function in 
each box with a parabola in the direction of the quarter chord of 
each panel, divided by r^^, where r^^ = + 2 .^ , y^ and z^ being 

relative y and z distances of the midpoints of the quarter chords 
(sending points) and the three qiiarter chord (receiving points). 
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The force on the panel due to the pressure distribution is assumed to 

be equal to the force on the doublet line segment. The accuracy of 

the method was remarkable for aeroelastic purposes. However for 

large frequencies (reduced freq'uencies k > l) it was found that the 

method was not only not accurate , but also did not converge to the 

results obtained by piston theory. It was found that by approximating 

the kernel in each box with a quartic in the spanwise direction of the 

2 

quarter chord point of each panel divided by r^^ , the resiilts were 
tremendously improved for planar lifting surfaces for high fre- 
quencies . 

Using the same notation as Hodden, Giesing and Kalman [ 77 ] we 
may rewrite equation 27b of Ref . 77 Si'S , 

P (q) [Kj^exp{-io)(x-n tan \ s)/v} - K-lq] 

(y-q)^ (y - q)^ 

- A^q + B^q + _ 2 

(y - q) 

Assuming that P (q) is evaluated at q = + ae , + be and 0, where e is 
the semi -span of each panel we have, 

A = 2 — -^(P(ae)+P(-ae)-2P(0))- -^(P(be)+P(-be)-2P(0) ) 

-t )e a b 

B. = 2A^ + J- p r i(P(ae)-P(-ae))- i(P(be)-P(-be)) 
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^(P(ae)-P(-ae))- -4(P(be)-P(-be)) 

ar J 

= P(0) 

*a^ and were chosen to be 1*0 and 0.5 respectively- For 
Aq = ®o ^ formulas for and reduce to those given 

by Albano and Rodden- The above approximation is very good for high 
frequencies, especially for swept wings (Xs ^ O) • At high frequencies, 
for swept wings, P (t]) varies rapidly for different values of How- 
ever extensive numerical experiments have to be carried out determine 
the optimm values of a, b and the receiving and sending points in 
various frequency domains - 

In spite of the improvement in the results for moderately high 
frequencies, the results did not converge to the piston theoiy limit 
for certain mode shapes - For convergence the distance between the 
sending and receiving point on each panel must vanish as this condition 
is known to hold in the case of piston theory. This would require, 
theoritically , an infinite number of boxes in the chordwise direction. 
Also for a proper spanwise pressure distribution, which critically 
effects the aerodynamic loads a large nxmiber of spanwise boxes are 
reqiiired at all frequencies. This leads to an excessively large 
computer storage requirements and thus the method is not too attrac- 
tive for high frequencies . 


B, = -B^y 


-2 


+ 




2(b^-a^)e 
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Laschka's kernel function method: In this method the pressure 

distribution is taken in the form 

40j,(x,y) . ^E^a^(y)hj(x) 

Following Laschka, the terms in the kernel function which are singular 
are separated from the part that is regular. The singular terms are 
integrated analytically while the regular terms are integrated 
numerically. The integral equation is reduced to a linear system and 
the values of the unknown coefficients a^(y) are obtained at different 
spanwise and chordwise stations • For the chordwise integration, the 
integrations points and the collocation points may be chosen inde- 
pendentally. For hi^ frequencies one may choose fewer collocation 
points eind a large number of integration points, in order to converge 
to piston theory. Thus the kernel function method is more efficient 
than the doublet lattice method. 

Control stir faces may be handled by both methods, without taking 
into account the singularities at their leading edges. However, the 
computational time for high frequencies is quite high when compared 
to the corresponding time at low frequencies, in both methods. 
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APPENDIX D 


PADE APPROXIMANTS AND INTEGRAL EQUATIONS 
Pad^ Approximants 

The [N,M] Pade approximant of a matrix function F(z) is defined hy 

F(z) = P(z)r"^(z) 


where 

P(z) is a matrix polynomial of degree M and 
r(z) is a matrix polynomial of degree N 

If in particular F(z) is a power series defined by^, 

L ^ 

F(z) = E F. z^ and if 
i=0 ^ 

M . i A 

P(z) = E P z^, R(z) = E R z , Rj^= I 
i=0 i=0 

we have the following equations for P^ and R^ 


N-1 


E 

J=0 



R. + F 


k-N 



k = 0, 1, 


M 


N-1 

E F, .R. = -F, k = M+ l,M+2,M+5, --.L 
j=o 


(D-5) 
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For L > M + N, equations D-5 may be solved for the matrices R. and P, 

— j K 

successively. 

On the other hand, if F(z) is known for z = k = 0,1, 2, 5 ••• L 
we have 

M . N-1 . 

F(z, )zj^ = E P zj,^ - F(z ) E R ^ k = 0,1,2,5,... L 
^ ^ i=0 i=0 

(D.4) 

It is now assumed that P^ and Pj^ are known apriori. (if Pj^ is not 
known, we may increase value of M by 1 and assume that = [O])- 

From the L known values of z.^, we pick M-1 values, 

— M-1, and rewrite equations D.4 in part as 


M-1 

E 

i=l 



N 




- p - p 

0 M k. 


+ F(z 



kj^ = 1, 2, 5, M-1 


We now have M-1 unknowns P^ and M-1 equations and hence we may solve 
for the vmknowns P^ in terms of R^ as 


P. 

1 


N-1 

E 

J=0 


B. .R . + 
ij 0 0 


(D.6) 


Hence equations D*4 reduce to a system of L-Mfl equations for the N 
unknowns R., j = 0,1,2, ... N-1, which may determine by least squares 
provided the system of eqmtions are solvable • On, the other hand it is 
also possible to solve equation D-4 directly by least squares. These 
equations are also applicable for scalar quantities . Properties of 
Fade approximants are disciissed fiirther in [78]- 
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Properties of Integral Equations delated to Fade Approximants 

([79] and [8O]) 

In this section some of the basic definitions and properties of 

I 

integral equations of the form 

_b 

cp(x) = f(x,s) + s J‘ dy K(x,y,s)q>(y) (D.7) 

a 

where f(x,s) and K(x,y,s) are defined on the closed internal [a^b], 
are discussed. We can regard ^ 

b 

K = J* dy K(x,y,s)(-) (D.8) 

a 

in equation D -7 as an operator acting on the function cp. In operator 
form the integral equation is 

cp = f + sKcp (D- 9 ) 

The function K(x,y,s) is the kernel of the equation and in order 
to ensure a valid definition and a unique solution of D»7, various 
restrictions can be imposed on cp, f, and K. It is usual to assume 
cp and f are functions of bounded norm || * || in the Hilbert space 
of functions defined on the interval [a,b]; for example, 

J" lcp(x) l^dx <~ (D.IO) 

a 

A sequence of functions [f is bounded if 
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f II < B < « 

n " — 


(D.ll) 


Bounded operators . - The operator defined by D.8 Is bounded 
if there is a finite constant M, such that 

llKhll <Ml|h|| (D.12) 

for all functions h of finite norm. This is a different condition from 
the boundedness of the function K (x,y,s) given by |K(x,y,s) | < c < » . 
The norm of a bounded operator K is the least value of M satisfying 
D.12. The operator is unbomded if no finite M can be fo\ind such that 
(D.12) holds for all h. 

Kernels of finite rank.- A kernel K is of finite rank n 
— n 

if it is expressible in the form 

K^(x,y,s) = s a^(x)P^(y) (D-15) 

i=l 

2 

where and are linearly independent sets of fxmctions in L • 

If the kernel K(x,y,s) is a kernel of finite rank and and are 
rationed, functions or polynomials of s eqiiation (D-7) may solved 
exactly. For simplicity it is assvaned that f(x,s) is independent of 

s . 

The kernel function method ; - The scalar product of functions 
2 

f emd g in L is 
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(D.14) 


b 

(f ^g) = J‘ (3y f (y)g(y) 

a 

Substituting (D-15) in (D-7) gives, 


n 

cp(x) = f(x) + s E (cp,PJa. (x) 


(D.I5) 


The coefficients (cp,P^) are complex numbers and einy solution of (D-7) 
may be of the form. 


n 

q) = f (x) + 2 F a 

d=l ^ ^ 


■where the coefficients and to be determined. 

Substituting (D-16) in (D.15), §. are determined by 


n 


i=l 




Let A(s) be the determinant of the coefficients of and the 

J 

minor associated with i,J^^ element be • 


Then 


.-1 


n 


and 




n n 

cp = f(x) + E H A s(f,S ) 

J=1 i=l 


,-l. 


n n 


Hence, cp = A' [Af(x)+ E E A . s (f (x) ,P. (x) )g ] 

d=l i=l '' 1 J 


(d.i6) 


As a function of the parameter s, this a rational function, the numera- 
tor degree and the denominator degree being dependent on and . 
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Compact operators . - Let T be a linear operator that trans- 

2 2 

forms £iny function h in L into another function Th of L ; in parti- 
cular T could he of the form given hy equation D.8. Every linear 
operator of this type is continuous in the sense that it transfers a 
sequence {h^^} converging in the mean to h into a sequence {fh^j} 
converging to Th. 

Compactness or complete continuity of a linear operator is a 
stronger condition than continuity. Let {f^} he any hounded infinite 

P 2 

sequence in L , so that CTfjj] is another infinite sequence in L . 
Then if for all hounded contains a subsequence converging 

to some function h as n so that 

II Tf^ - hll =0 


then T is completely continuous. Reiz and St. -Nagy [80 ] give necessary 
and sorfficient conditions for complete continuity . 

Completely continuous operators can he uniformly approximated hy 
operators with kernels of finite rank. Also, if a kernel function 
can he approximated with kernels of finite rank, then the operator 
is completely continuous . In unsteady aerodynamics the kernels are 
implicitly approximated hy kernels of finite rank and resulting solu- 
tions are known to converge to the exact press\are distributions vini- 
formly, for all frequencies . 

Thus there exists an infinite sequence of kernels ... 

of rank 1 , 2 , 3 , • • • such that 
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II (K - K^)h|| < e^llhll 


with _> 0 as n ^ “ for all h € uniformly for all s . 


m 



are rational functions or polynomials in the variables s. Hence it 
may be assumed that as s ® is of O(s^) for all n. Then from 
equations D-l6 we obtain a sequence of solutions, which are rational 
functions of s . The degree of the numerator would be W-E that of 
“the d6nom.ina.tor N \diere N a.s n . 

Neumann series solutions . - 

Another method of solving (D-T) is to substitute the series 

00 

cp(x) = f(x) + E s“fj^(x,s) (D.I7) 

n=l 

and equate the coefficients of s on each side . This gives 

b 

^n+1^^'®^ = J dy K(x,y,s)fn(y,s) 
a 

n = 1,2,5, ... 

The series (D-l?) is known as a series and may be written in 

operator notation as , 

cp(x) = f(x) + Kf(x) + K^f(x) ... 


When K is uniformly bounded in the norm, then the series converges in 
the norm if. 
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-1 


Isl < 11 k 11 

When s does not satisfy this condition Chisholm [79] has shown that 
Fade approximants of the series may be formed which will give approxi- 
mations to the solution for all s • 

From Chisholm's theorem it follows that: If the kernel 

function K(x, y, s)of the integral equation D-7 can be approximated 
by a sequence of kernels of finite rahk^ then there exists a sequence of 
solutions which are rational functions in the 's' domain, the degree of the 
numerator being N + J and the degree of the denominator being N, which 
converges uniformly to the exact solution as N > where J depends 
on the limiting behavior of K(x,y,s) and f(x,s) as s -v “ • 

Thus if the behavior of the pressure distributions in unsteady 
aerodynamics are known for high frequencies and in steady flow and if 
it is possible to show that the associated kernel functions are 
completely continuous for all s, then it is possible, in principle, 
to construct Fade approximants for the aerodynamic loads from the 
numerically obtained solutions • 
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1 refers to pltinging, 2 to pitching about mid chord , and 3 to the flap mode. 



TABLE 2. 

Coefficients of Rational F\mction Approximations to Non-dimensional Aerodynamic Load Coefficients 
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TABLE 5. 

Coefficients of Rational Function Approximations to Non-dimensional Aerodynamic Load Coefficients 

for a Two-dimensional Airfoil, M = 0.5 
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TABLE 4. 

Coefficients of Rational Function Approximations to Non-dimensional Aerodynamic Load Coefficients 
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Coefficients of Rational Function Approximations to Non-dimensional Aerodynamic Load Coefficients 

for a Two-dimensional Airfoil, M = 0.7 
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TABLE 6. 

Coefficients of RationaJ. Function Approximations to Non-dimensional Aerodynamic Load Coefficients 

for a Two-dimensional Airfoil, M = 1.5 
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Coefficients of Rational Function Approximations to Non-dimensional Aerodynamic Load Coefficients 

for a Two-dimensional Airfoil, M = 2.0 
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Rational Function Approximations of Generalized Aerodynamic Loads in Plunging (l). 
Pitching ( 2 ), Chordwise Bending (x/r) , (5) and, Spanwise Bending (y/rj (ii-) 
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Indicial Lift in Plvinging for circular and Elliptic Wings 
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Fig- 4. Error in F(k) (True - Approximate) 




Fig- 5* Error in G(k) (True - Approximate). 
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distance traveled in semi-chords ^ 
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Fig- 11- Indicial lift for an impulsive plunging displacement^ AR 
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distance traveled in semi-chords 

Fig- 17- Indicial partial moment in plunging. 
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•asymptotes, low aspect 



Fig. 26. Variation of C with aspect ratio for rectangular and elliptic wings 





Pi O 


5 ^ 


o 

Vi 

-P P 

V cd 


O Pi 
•H *H 
■d rH 
P rH 


t^ 

OJ 

bO 

•H 

pH 


160 


asymptotes 
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Fig. 28. Indicial lift due to impulsive 





distance traveled in semi-chords T - 

Indicial lift due to impulsive plunging^ rectangular planform, AR 
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Fig. 52. Influence of finite aspect ratio and sweep on the indicial function for impulsive 
plunging at M = 0-7- 
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due to impulsive displacement in the first bending mode, 
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Fig. Indicial function for due to impulsive flap rotation. 
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Fig. 4i+. Indicial function for q 
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Control Laws For Flutter 


Suppression 


Fig. 47- Block diagram for control law synthesis procedure. 
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